Changeset 1a8fbee in flexpart.git for src/verttransform_ecmwf.f90


Ignore:
Timestamp:
Jun 17, 2018, 5:14:02 PM (6 years ago)
Author:
pesei <petra seibert at univie ac at>
Branches:
univie
Children:
8b3d324
Parents:
77778f8
Message:

vertransform: fix #140 (ref position for grid), verttransform_nest: fix bug with dimensions, both: smaller changes, com_mod: correct comment to nmixz

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/verttransform_ecmwf.f90

    r47f96e5 r1a8fbee  
    1 !**********************************************************************
     1! **********************************************************************
    22! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    … …  
    1818! You should have received a copy of the GNU General Public License   *
    1919! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
    20 !**********************************************************************
     20! **********************************************************************
    2121
    2222subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
    … …  
    3636!     Update: 16 January 1998                                                *
    3737!                                                                            *
     38!                                                                            *
     39!*****************************************************************************
     40!  CHANGES                                                                   *
    3841!     Major update: 17 February 1999                                         *
    3942!     by G. Wotawa                                                           *
    … …  
    4346!       procedure                                                            *
    4447!                                                                            *
    45 !*****************************************************************************
    46 !  Changes, Bernd C. Krueger, Feb. 2001:
     48!  Bernd C. Krueger, Feb. 2001:
    4749!   Variables tth and qvh (on eta coordinates) from common block
    48 !
    49 ! Sabine Eckhardt, March 2007
    50 ! added the variable cloud for use with scavenging - descr. in com_mod
    51 !
    52 ! Unified ECMWF and GFS builds
     50!                                                                            *
     51! Sabine Eckhardt, March 2007:
     52!   added the variable cloud for use with scavenging - descr. in com_mod
     53!                                                                            *
     54! Who? When?                                                                 *
     55!   Unified ECMWF and GFS builds
    5356! Marian Harustak, 12.5.2017
    5457!     - Renamed from verttransform to verttransform_ecmwf
    55 !*****************************************************************************
    56 ! Date: 2017-05-30 modification of a bug in ew. Don Morton (CTBTO project)   *
     58!                                                                            *
     59! Don Morton, 2017-05-30:
     60!   modification of a bug in ew. Don Morton (CTBTO project)                  *
     61!                                                                            *
     62!  undocumented modifications by NILU for v10                                *
     63!                                                                            *
     64!  Petra Seibert, 2018-06-13:                                                *
     65!   - fix bug of ticket:140 (find reference position for vertical grid)      *
     66!   - put back SAVE attribute for INIT, just to be safe                      *
     67!   - minor changes, most of them just cosmetics                             *
     68!   for details see changelog.txt                                            *
     69!                                                                            *
     70! ****************************************************************************
     71
    5772!*****************************************************************************
    5873!                                                                            *
    … …  
    7388  use com_mod
    7489  use cmapf_mod, only: cc2gll
    75 !  use mpi_mod
    7690
    7791  implicit none
    … …  
    8296  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
    8397  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
     98!>   for reference profile (PS) 
     99  real ::  tvoldref, poldref, pintref, psmean, psstd
     100  integer :: ixyref(2)
     101  integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you
     102!!   want old method of searching reference location for the vertical grid
     103!!   1 for new method (options for other methods 2,... in the future)
     104
     105  real,dimension(0:nymax-1) :: cosf
    84106  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
    85   real,dimension(0:nymax-1) :: cosf
     107!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
    86108
    87109  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
    88110
    89   integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
     111  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixref,jyref,kz_inv
    90112  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
    91113  real :: ew,dz1,dz2,dz
    … …  
    96118  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
    97119
    98   logical :: init = .true.
    99   logical :: init_w = .false.
    100   logical :: init_r = .false.
     120  logical, save :: init = .true. ! PS 2018-06-13: add back save attribute
    101121
    102122
    … …  
    120140! z levels in meter. The heights are the heights of model levels, where  *
    121141! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
    122 ! the vertical resolution in the z system is doubled. As reference point,*
     142! the vertical rESOlution in the z system is doubled. As reference point,*
    123143! the lower left corner of the grid is used.                             *
    124144! Unlike in the eta system, no difference between heights for u,v and    *
    … …  
    126146!*************************************************************************
    127147
    128 
    129 !eso measure CPU time
     148!ESO measure CPU time
    130149!  call mpif_mtime('verttransform',0)
    131150
    132151  if (init) then
    133152
    134 
    135     if (init_r) then
    136 
    137         open(333,file='heights.txt', &
    138           form='formatted')
    139         do kz=1,nuvz
    140             read(333,*) height(kz)
    141         end do
    142         close(333)
    143         write(*,*) 'height read'
    144     else
    145 
    146 
    147 ! Search for a point with high surface pressure (i.e. not above significant topography)
    148 ! Then, use this point to construct a reference z profile, to be used at all times
    149 !*****************************************************************************
    150 
    151     do jy=0,nymin1
    152       do ix=0,nxmin1
    153         if (ps(ix,jy,1,n).gt.100000.) then
    154           ixm=ix
    155           jym=jy
    156           goto 3
    157         endif
    158       end do
     153!! Search for a point with high surface pressure (i.e. not above significant
     154!!  topography) Then, use this point to construct a reference z profile,
     155!!  to be used at all times
     156! *****************************************************************************
     157
     158    if (ioldref .eq. 0) then ! old reference grid point
     159       do jy=0,nymin1
     160         do ix=0,nxmin1
     161           if (ps(ix,jy,1,n).gt.1000.e2) then
     162             ixref=ix
     163             jyref=jy
     164             goto 3
     165           endif
     166         end do
     167       end do
     1683     continue
     169
     170      print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
     171    else ! new reference grid point
     172!     PS: the old version fails if the pressure is <=1000 hPa in the whole
     173!     domain. Let us find a good replacement, not just a quick fix.
     174!     Search point near to mean pressure after excluding mountains
     175
     176      psmean = sum( ps(:,:,1,n) ) / (nx*ny)
     177      print*,'height: fg psmean',psmean
     178      psstd = sqrt(sum( (ps(:,:,1,n) - psmean)**2 ) / (nx*ny))
     179
     180!>    new psmean using only points within $\plusminus\sigma$ 
     181!      psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / &
     182!        count(abs(ps(:,:,1,n)-psmean) < psstd)
     183
     184!>    new psmean using only points with $p\gt \overbar{p}+\sigma_p$ 
     185!>    (reject mountains, accept valleys)
     186      psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / &
     187        count(ps(:,:,1,n) > psmean - psstd)
     188      print*,'height: std, new psmean',psstd,psmean
     189      ixyref = minloc( abs( ps(:,:,1,n) - psmean ) )
     190      ixref = ixyref(1)
     191      jyref = ixyref(2)
     192      print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
     193    endif 
     194
     195    tvoldref=tt2(ixref,jyref,1,n)* &
     196      ( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n))
     197    poldref=ps(ixref,jyref,1,n)
     198    height(1)=0.
     199
     200    do kz=2,nuvz
     201
     202      pintref = akz(kz)+bkz(kz)*ps(ixref,jyref,1,n)
     203      tv = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n))
     204
     205      if (abs(tv(ixref,jyref)-tvoldref) .gt. 0.2) then
     206        height(kz) = height(kz-1) +      &
     207          const*log( poldref/pintref ) * &
     208          ( tv(ixref,jyref) - tvoldref ) / log( tv(ixref,jyref) / tvoldref )
     209      else
     210        height(kz) = height(kz-1) + &
     211          const*log( poldref/pintref ) * tv(ixref,jyref)
     212      endif
     213
     214      tvoldref = tv(ixref,jyref)
     215      poldref = pintref
     216      print*,'height=',kz,height(kz),tvoldref,poldref
     217
    159218    end do
    160 3   continue
    161 
    162 
    163     tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    164          ps(ixm,jym,1,n))
    165     pold(ixm,jym)=ps(ixm,jym,1,n)
    166     height(1)=0.
    167 
    168     do kz=2,nuvz
    169       pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
    170       tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    171 
    172       if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
    173         height(kz)= &
    174              height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
    175              (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
    176       else
    177         height(kz)=height(kz-1)+ &
    178              const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
    179       endif
    180 
    181       tvold(ixm,jym)=tv(ixm,jym)
    182       pold(ixm,jym)=pint(ixm,jym)
    183     end do
    184 
    185     if (init_w) then
    186         open(333,file='heights.txt', &
    187           form='formatted')
    188         do kz=1,nuvz
    189               write(333,*) height(kz)
    190         end do
    191         close(333)
    192     endif
    193 
    194     endif ! init
     219
    195220
    196221! Determine highest levels that can be within PBL
    … …  
    207232! Do not repeat initialization of the Cartesian z grid
    208233!*****************************************************
    209 
    210234    init=.false.
    211235
    212 !    dbg_height = height
    213 
    214   endif
     236  endif ! init block (vertical grid construction)
    215237
    216238
    … …  
    218240!*************************
    219241
    220 
    221   do jy=0,nymin1
    222     do ix=0,nxmin1
    223       tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
    224            ps(ix,jy,1,n))
    225     enddo
    226   enddo
     242  tvold(:,:)=tt2(:,:,1,n) * (1.+0.378*ew( td2(:,:,1,n) ) / ps(:,:,1,n))
    227243  pold=ps(:,:,1,n)
    228244  uvzlev(:,:,1)=0.
    … …  
    235251
    236252  do kz=2,nuvz
    237     pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
    238     tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
    239     rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
    240 
    241     where (abs(tv-tvold).gt.0.2)
    242       uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
    243            (tv-tvold)/log(tv/tvold)
     253    pint(:,:)=akz(kz)+bkz(kz)*ps(:,:,1,n)
     254    tv(:,:)=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
     255    rhoh(:,:,kz)=pint(:,:)/(r_air*tv(:,:))
     256
     257    where (abs(tv(:,:)-tvold(:,:)).gt.0.2)
     258      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))* &
     259           (tv(:,:)-tvold(:,:))/log(tv(:,:)/tvold(:,:))
    244260    elsewhere
    245       uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
     261      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))*tv(:,:)
    246262    endwhere
    247263
    248     tvold=tv
    249     pold=pint
     264    tvold(:,:)=tv(:,:)
     265    pold(:,:)=pint(:,:)
    250266  end do
    251267
    … …  
    260276
    261277  pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
    262        ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
    263        (aknew(1)+bknew(1)*ps(:,:,1,n)))
     278    ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
     279     (aknew(1)+bknew(1)*ps(:,:,1,n)))
    264280  do kz=2,nz-1
    265281    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
    266          ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
    267          (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
     282      ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
     283       (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
    268284  end do
    269285  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
    270286       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
    271        (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
    272 
    273 ! Levels, where u,v,t and q are given
    274 !************************************
    275 
     287        (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
     288
     289! Levels where u,v,t and q are given
     290!***********************************
    276291
    277292  uu(:,:,1,n)=uuh(:,:,1)
    … …  
    279294  tt(:,:,1,n)=tth(:,:,1,n)
    280295  qv(:,:,1,n)=qvh(:,:,1,n)
    281 !hg adding the cloud water
     296!HG adding the cloud water
    282297  if (readclouds) then
    283298    clwc(:,:,1,n)=clwch(:,:,1,n)
    284299    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
    285300  end if
    286 !hg
     301!HG
    287302  pv(:,:,1,n)=pvh(:,:,1)
    288303  rho(:,:,1,n)=rhoh(:,:,1)
    … …  
    292307  tt(:,:,nz,n)=tth(:,:,nuvz,n)
    293308  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
    294 !hg adding the cloud water
     309!HG adding the cloud water
    295310  if (readclouds) then
    296311    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
    297     if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
     312    if (.not. sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
    298313  end if
    299 !hg
     314!HG
    300315  pv(:,:,nz,n)=pvh(:,:,nuvz)
    301316  rho(:,:,nz,n)=rhoh(:,:,nuvz)
    302317
    303 
    304318  kmin=2
    305319  idx=kmin
    306   do iz=2,nz-1
     320  iz_loop: do iz=2,nz-1
    307321    do jy=0,nymin1
    308322      do ix=0,nxmin1
    309         if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
     323        if (height(iz).gt.uvzlev(ix,jy,nuvz)) then
     324       
    310325          uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    311326          vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    312327          tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
    313328          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
    314 !hg adding the cloud water
     329!HG adding the cloud water
    315330          if (readclouds) then
    316331            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
    317             if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     332            if (.not. sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
    318333          end if
    319 !hg
     334!HG
    320335          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
    321336          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
     337
    322338        else
    323           innuvz: do kz=idx(ix,jy),nuvz
    324             if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    325                  (height(iz).le.uvzlev(ix,jy,kz))) then
     339
     340          kz_loop: do kz=idx(ix,jy),nuvz
     341            if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
     342                .and. height(iz).le.uvzlev(ix,jy,kz)) then
    326343              idx(ix,jy)=kz
    327               exit innuvz
     344              exit kz_loop
    328345            endif
    329           enddo innuvz
     346          enddo kz_loop
     347
    330348        endif
    331349      enddo
    … …  
    333351    do jy=0,nymin1
    334352      do ix=0,nxmin1
    335         if(height(iz).le.uvzlev(ix,jy,nuvz)) then
     353        if (height(iz).le.uvzlev(ix,jy,nuvz)) then
     354       
    336355          kz=idx(ix,jy)
    337356          dz1=height(iz)-uvzlev(ix,jy,kz-1)
    … …  
    344363          qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
    345364               +qvh(ix,jy,kz,n)*dz1)/dz
    346 !hg adding the cloud water
     365!HG adding the cloud water
    347366          if (readclouds) then
    348367            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
    349             if (.not.sumclouds) &
    350                  &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
     368            if (.not. sumclouds) ciwc(ix,jy,iz,n)= &
     369              (ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
    351370          end if
    352 !hg
     371!HG
    353372          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    354373          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
     374
    355375        endif
    356376      enddo
    357377    enddo
    358   enddo
    359 
    360 
    361 ! Levels, where w is given
     378  enddo iz_loop
     379
     380
     381! Levels where w is given
    362382!*************************
    363383
    … …  
    366386  kmin=2
    367387  idx=kmin
    368   do iz=2,nz
     388  iz_loop2: do iz=2,nz
    369389    do jy=0,nymin1
    370390      do ix=0,nxmin1
    371         inn:        do kz=idx(ix,jy),nwz
    372           if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
    373               height(iz).le.wzlev(ix,jy,kz)) then
     391        kz_loop2: do kz=idx(ix,jy),nwz
     392          if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
     393            .and. height(iz).le.wzlev(ix,jy,kz)) then
    374394            idx(ix,jy)=kz
    375             exit inn
     395            exit kz_loop2
    376396          endif
    377         enddo inn
     397        enddo kz_loop2
    378398      enddo
    379399    enddo
    … …  
    385405        dz=dz1+dz2
    386406        ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
    387              +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
    388       enddo
    389     enddo
    390   enddo
     407                      + wwh(ix,jy,kz)  *pinmconv(ix,jy,kz)*dz1)/dz
     408      enddo
     409    enddo
     410  enddo iz_loop2
    391411
    392412! Compute density gradients at intermediate levels
    393413!*************************************************
    394414
    395   drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
    396        (height(2)-height(1))
     415  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/(height(2)-height(1))
    397416  do kz=2,nz-1
    398417    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
    399          (height(kz+1)-height(kz-1))
     418      (height(kz+1)-height(kz-1))
    400419  end do
    401420  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
    402 
    403 !    end do
    404 !  end do
    405421
    406422
    … …  
    411427
    412428  do jy=1,ny-2
    413     cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
     429    cosf(jy) = 1. / cos( ( real(jy)*dy + ylat0 )*pi180 )
    414430  enddo
    415431
    416432  kmin=2
    417433  idx=kmin
    418   do iz=2,nz-1
     434  iz_loop3: do iz=2,nz-1
    419435    do jy=1,ny-2
    420436      do ix=1,nx-2
    421437
    422         inneta: do kz=idx(ix,jy),nz
    423           if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    424                (height(iz).le.uvzlev(ix,jy,kz))) then
     438        kz_loop3: do kz=idx(ix,jy),nz
     439          if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
     440            .and. height(iz).le.uvzlev(ix,jy,kz)) then
    425441            idx(ix,jy)=kz
    426             exit inneta
     442            exit kz_loop3
    427443          endif
    428         enddo inneta
     444        enddo kz_loop3
    429445      enddo
    430446    enddo
    … …  
    442458
    443459        dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
    444         dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
     460        dzdx2=(uvzlev(ixp,jy,kz)  -uvzlev(ix1,jy,kz)  )/2.
    445461        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
    446462
    447463        dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
    448         dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
     464        dzdy2=(uvzlev(ix,jyp,kz)  -uvzlev(ix,jy1,kz)  )/2.
    449465        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    450466
    451         ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst)
    452 
    453       end do
    454 
    455     end do
    456   end do
     467        ww(ix,jy,iz,n)=ww(ix,jy,iz,n) +( dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy) &
     468          + dzdy*vv(ix,jy,iz,n)*dyconst)
     469
     470      enddo
     471    enddo
     472
     473  enddo iz_loop3
    457474
    458475! If north pole is in the domain, calculate wind velocities in polar
    … …  
    461478
    462479  if (nglobal) then
     480   
    463481    do iz=1,nz
    464482      do jy=int(switchnorthg)-2,nymin1
    … …  
    467485          xlon=xlon0+real(ix)*dx
    468486          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    469                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    470                vvpol(ix,jy,iz,n))
    471         end do
    472       end do
    473     end do
     487               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), vvpol(ix,jy,iz,n))
     488        enddo
     489      enddo
     490    enddo
    474491
    475492
    476493    do iz=1,nz
    477494
    478 ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
     495!   CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
    479496!
    480 !   AMSnauffer Nov 18 2004 Added check for case vv=0
    481 !
     497!   AMS nauffer Nov 18 2004 Added check for case vv=0
     498
    482499      xlon=xlon0+real(nx/2-1)*dx
    483500      xlonr=xlon*pi/180.
    484       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
    485            vv(nx/2-1,nymin1,iz,n)**2)
     501      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
    486502      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
    487         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
    488              vv(nx/2-1,nymin1,iz,n))-xlonr
     503        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
    489504      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    490         ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    491              vv(nx/2-1,nymin1,iz,n))-xlonr
     505        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
    492506      else
    493         ddpol=pi/2-xlonr
     507        ddpol=pi/2.-xlonr
    494508      endif
    495       if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
    496       if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
     509      if (ddpol.lt.0.) ddpol=2.0*pi+ddpol
     510      if (ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
    497511
    498512! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
    … …  
    502516      uuaux=-ffpol*sin(xlonr+ddpol)
    503517      vvaux=-ffpol*cos(xlonr+ddpol)
    504       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    505            vvpolaux)
     518      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, vvpolaux)
    506519
    507520      jy=nymin1
    … …  
    509522        uupol(ix,jy,iz,n)=uupolaux
    510523        vvpol(ix,jy,iz,n)=vvpolaux
    511       end do
    512     end do
    513 
    514 
    515 ! Fix: Set W at pole to the zonally averaged W of the next equator-
    516 ! ward parallel of latitude
     524      enddo
     525    enddo
     526
     527! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next
     528! equator-ward parallel
    517529
    518530    do iz=1,nz
    … …  
    521533      do ix=0,nxmin1
    522534        wdummy=wdummy+ww(ix,jy,iz,n)
    523       end do
     535      enddo
    524536      wdummy=wdummy/real(nx)
    525537      jy=nymin1
    526538      do ix=0,nxmin1
    527539        ww(ix,jy,iz,n)=wdummy
    528       end do
    529     end do
     540      enddo
     541    enddo
    530542
    531543  endif
    … …  
    537549
    538550  if (sglobal) then
     551 
    539552    do iz=1,nz
    540553      do jy=0,int(switchsouthg)+3
    … …  
    543556          xlon=xlon0+real(ix)*dx
    544557          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    545                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    546                vvpol(ix,jy,iz,n))
    547         end do
    548       end do
    549     end do
     558               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     559        enddo
     560      enddo
     561    enddo
    550562
    551563    do iz=1,nz
    552564
    553 ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
     565!   CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
    554566!
    555 !   AMSnauffer Nov 18 2004 Added check for case vv=0
    556 !
     567!   AMS nauffer Nov 18 2004 Added check for case vv=0
     568
    557569      xlon=xlon0+real(nx/2-1)*dx
    558570      xlonr=xlon*pi/180.
    559       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
    560            vv(nx/2-1,0,iz,n)**2)
     571      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
    561572      if (vv(nx/2-1,0,iz,n).lt.0.) then
    562         ddpol=atan(uu(nx/2-1,0,iz,n)/ &
    563              vv(nx/2-1,0,iz,n))+xlonr
     573        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    564574      else if (vv(nx/2-1,0,iz,n).gt.0.) then
    565         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
    566              vv(nx/2-1,0,iz,n))+xlonr
     575        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    567576      else
    568         ddpol=pi/2-xlonr
     577        ddpol=pi/2.-xlonr
    569578      endif
    570       if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
    571       if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
     579      if (ddpol.lt.0.) ddpol=2.0*pi+ddpol
     580      if (ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
    572581
    573582! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
    … …  
    577586      uuaux=+ffpol*sin(xlonr-ddpol)
    578587      vvaux=-ffpol*cos(xlonr-ddpol)
    579       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    580            vvpolaux)
     588      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    581589
    582590      jy=0
    … …  
    584592        uupol(ix,jy,iz,n)=uupolaux
    585593        vvpol(ix,jy,iz,n)=vvpolaux
    586       end do
    587     end do
    588 
    589 
    590 ! Fix: Set W at pole to the zonally averaged W of the next equator-
    591 ! ward parallel of latitude
     594      enddo
     595    enddo
     596
     597! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next
     598! equator-ward parallel
    592599
    593600    do iz=1,nz
    … …  
    596603      do ix=0,nxmin1
    597604        wdummy=wdummy+ww(ix,jy,iz,n)
    598       end do
     605      enddo
    599606      wdummy=wdummy/real(nx)
    600607      jy=0
    601608      do ix=0,nxmin1
    602609        ww(ix,jy,iz,n)=wdummy
    603       end do
    604     end do
     610      enddo
     611    enddo
    605612  endif
    606613
    607614
    608 !*********************************************************************************** 
    609   if (readclouds) then !HG METHOD
    610 ! The method is loops all grids vertically and constructs the 3D matrix for clouds
    611 ! Cloud top and cloud bottom gid cells are assigned as well as the total column
    612 ! cloud water. For precipitating grids, the type and whether it is in or below
     615!******************************************************************************
     616  if (readclouds) then ! HG METHOD
     617
     618! Loops over all grid cells vertically and construct the 3D matrix for clouds
     619! Cloud top and cloud bottom grid cells are assigned as well as the total column
     620! cloud water. For precipitating columns, the type and whether it is in or below
    613621! cloud scavenging are assigned with numbers 2-5 (following the old metod).
    614 ! Distinction is done for lsp and convp though they are treated the same in regards
    615 ! to scavenging. Also clouds that are not precipitating are defined which may be
    616 ! to include future cloud processing by non-precipitating-clouds.
    617 !***********************************************************************************
    618     write(*,*) 'Global ECMWF fields: using cloud water'
     622! A distinction is made between lsp and convp though they are treated the equally
     623! with regard to scavenging. Also, clouds that are not precipitating are defined which
     624! may be used in the future for cloud processing by non-precipitating-clouds.
     625!*******************************************************************************
     626
     627!PS    write(*,*) 'Global ECMWF fields: using cloud water'
    619628    clw(:,:,:,n)=0.0
    620 !    icloud_stats(:,:,:,n)=0.0
    621629    ctwc(:,:,n)=0.0
    622630    clouds(:,:,:,n)=0
    623631! If water/ice are read separately into clwc and ciwc, store sum in clwc
    624     if (.not.sumclouds) then
     632    if (.not. sumclouds) then
    625633      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
    626     end if
     634    endif
    627635    do jy=0,nymin1
    628636      do ix=0,nxmin1
    629637        lsp=lsprec(ix,jy,1,n)
    630638        convp=convprec(ix,jy,1,n)
    631         prec=lsp+convp
     639        prec=lsp+convp  ! Note PS: prec is not used currently
    632640!        tot_cloud_h=0
    633641! Find clouds in the vertical
     642!! Note PS: bad loop order.
    634643        do kz=1, nz-1 !go from top to bottom
    635644          if (clwc(ix,jy,kz,n).gt.0) then     
    636645! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
    637             clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
     646            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))* &
     647              (height(kz+1)-height(kz))
    638648!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
    639649           
    … …  
    651661
    652662! If Precipitation. Define removal type in the vertical
    653         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     663        if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
     664!! Note PS: such hardcoded limits would better be parameters defined in par_mod
    654665
    655666          do kz=nz,1,-1 !go Bottom up!
     667!! Note PS: bad loop order
    656668            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
    657669              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
    … …  
    662674                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
    663675              endif                                              ! convective or large scale
    664             elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
     676            elseif( clw(ix,jy,kz,n).le.0 .and. cloudh_min.ge.height(kz)) then ! is below cloud
    665677              if (lsp.ge.convp) then
    666678                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
    … …  
    671683
    672684            if (height(kz).ge. 19000) then                        ! set a max height for removal
     685!! Note PS: such hardcoded limits would better be parameters defined in par_mod
    673686              clouds(ix,jy,kz,n)=0
    674             endif !clw>0
    675           end do !nz
     687            endif ! clw>0
     688          enddo ! kz
    676689        endif ! precipitation
    677690      end do
    678691    end do
    679692
    680 ! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
     693! ESO: copy the relevant data to clw4 to reduce amount of communicated data for MPI
    681694!    ctwc(:,:,n) = icloud_stats(:,:,4,n)
    682695
    … …  
    684697  else       ! use old definitions
    685698!**************************************************************************
     699
    686700!   create a cloud and rainout/washout field, clouds occur where rh>80%
    687701!   total cloudheight is stored at level 0
    688     write(*,*) 'Global fields: using cloud water from Parameterization'
     702
     703!PS    write(*,*) 'Global fields: using cloud water from Parameterization'
    689704    do jy=0,nymin1
    690705      do ix=0,nxmin1
    691 ! OLD METHOD
    692706        rain_cloud_above(ix,jy)=0
    693707        lsp=lsprec(ix,jy,1,n)
    … …  
    695709        cloudsh(ix,jy,n)=0
    696710        do kz_inv=1,nz-1
     711!! Note PS: bad loop order.
    697712          kz=nz-kz_inv+1
    698713          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    699714          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    700715          clouds(ix,jy,kz,n)=0
    701           if (rh.gt.0.8) then ! in cloud
    702             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     716
     717          if (rh .gt. 0.8) then ! in cloud
     718!! Note PS: such hardcoded limits would better be parameters defined in par_mod
     719
     720            if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
     721!! Note PS: such hardcoded limits would better be parameters defined in par_mod
    703722              rain_cloud_above(ix,jy)=1
    704723              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    … …  
    712731              clouds(ix,jy,kz,n)=1 ! cloud
    713732            endif
     733
    714734          else ! no cloud
     735
    715736            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
    716737              if (lsp.ge.convp) then
    … …  
    720741              endif
    721742            endif
     743
    722744          endif
    723745        end do
    724 !END OLD METHOD
     746
    725747      end do
    726748    end do
    727   endif !readclouds
     749  endif ! END OLD METHOD
    728750
    729751
    … …  
    869891!            icloudthck(ix,jy,n)=icmv
    870892!          endif
    871 !hg__________________________________
     893!HG__________________________________
    872894!           rcw(ix,jy)=2E-7*prec**0.36
    873895!           rpc(ix,jy)=prec
    874 !hg end______________________________
    875 
    876 !      endif !hg read clouds
    877 
    878 
    879 
    880 
    881 !eso measure CPU time
     896!HG end______________________________
     897
     898!      endif !HG read clouds
     899
     900
     901
     902
     903!ESO measure CPU time
    882904!  call mpif_mtime('verttransform',1)
    883905
    884 !eso print out the same measure as in Leo's routine
     906!ESO print out the same measure as in Leo's routine
    885907    ! write(*,*) 'verttransform: ', &
    886908    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG