Changeset 1a8fbee in flexpart.git for src


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

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/changelog.txt

    r77778f8 r1a8fbee  
    22Created by Petra Seibert, 8 June 2018 to document and explain in detail why I changed what.
    33
    4 makefile, 2018-06-08
    5 ========
     4* makefile, 2018-06-08
     5  ========
    66
    771) Add a GPL3+ License statement, author statement, version information. Prudent to do that, not harmful.
    … …  
    24243. What is the purpose of -Warray-bounds ? Is it relevant for Fortran? I think if there is a compile-time array-bound violation this would be a hard error?
    2525
    26 FLEXPART.f90
    27 ========
     26* FLEXPART.f90
     27  ========
    2828
    2929Update version string!
    3030
    31 readcommand.f90
    32 ===========
     31* readcommand.f90
     32  ===========
    3333
    3434Correct misleading error message (replace "open" by "write to")
    … …  
    3636Let STOP say in which subroutine we are stopping
    3737Some minor changes
     38
     39* verttransform_ecmwf.f90 PS 2018-06-13
     40  ===================
     41
     421) Remove some of commented out testing stuff
     432) Fix indenting in the init if block
     443) Code layout improvement
     454) Change name of loops to represent the index
     465) Fix ticket:140 Introduce a new way of establishing the reference position
     47   for the vertical grid. Also correct a minor bug in the calculation
     48   of height (array assignment where it was not intended)
     496) Add back the SAVE attribute to INIT, just to be sure
     507) Bring some more code under the IF (INIT) block
     518) make some annotations
     52
     53* com_mod.f90  PS 2018-06-13
     54  =======
     55 
     56  correct the comment for nmixz
     57
     58* verttransform_nest.f90 PS 2018-06-17
     59  ===================
     601) insert proper boundaries for implied loops in array expressions
     61  (fixes bug annoted by ESO in 2016)
     622) Code layout improvement
     633) Change name of loops to represent the index
     644) make some annotations

  • src/com_mod.f90

    rc2bd55e r1a8fbee  
    302302  ! nxfield                 same as nx for limited area fields,
    303303  !                    but for global fields nx=nxfield+1
    304   ! nmixz                   number of levels up to maximum PBL height (3500 m)
     304  ! nmixz                   number of levels up to maximum PBL height (hmixmax)
    305305
    306306  ! nuvz is used for u,v components

  • 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)), &

  • src/verttransform_nests.f90

    r72d3a5a r1a8fbee  
    3838!     Update: 16 January 1998                                                *
    3939!                                                                            *
     40!                                                                            *
     41!*****************************************************************************
     42!  CHANGES                                                                   *
    4043!     Major update: 17 February 1999                                         *
    4144!     by G. Wotawa                                                           *
    … …  
    4851!  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
    4952!   Variables tthn and qvhn (on eta coordinates) from common block
    50 !*****************************************************************************
    51 ! Sabine Eckhardt, March 2007
    52 ! add the variable cloud for use with scavenging - descr. in com_mod
    53 !*****************************************************************************
     53! Sabine Eckhardt, March 2007:
     54!   added the variable cloud for use with scavenging - descr. in com_mod
     55!                                                                            *
     56! Who? When?                                                                 *
     57!   Unified ECMWF and GFS builds
     58! Marian Harustak, 12.5.2017
     59!     - Renamed from verttransform to verttransform_ecmwf
     60!                                                                            *
    5461! ESO, 2016
    5562! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than
    56 !  the actual field dimensions
    57 !*****************************************************************************
    58 ! Date: 2017-05-30 modification of a bug in ew. Don Morton (CTBTO project)   *
    59 !*****************************************************************************
     63!  the actual field dimensions /fixed, PS 2018/
     64!                                                                            *
     65! Don Morton, 2017-05-30:
     66!   modification of a bug in ew. Don Morton (CTBTO project)                  *
     67!                                                                            *
     68!  undocumented modifications by NILU for v10                                *
     69!                                                                            *
     70!  Petra Seibert, 2018-06-13:                                                *
     71!   - insert proper boundaries for implied loops in array expressions        *s
     72!   - minor changes, most of them just cosmetics                             *
     73!   for details see changelog.txt                                            *
     74!                                                                            *
     75! ****************************************************************************
    6076!                                                                            *
    6177! Variables:                                                                 *
    … …  
    7591  implicit none
    7692
    77   real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn
     93  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: &
     94    uuhn,vvhn,pvhn
    7895  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
    7996
    8097  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
    8198  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
     99
     100  real,dimension(0:nymaxn-1) :: cosf
    82101  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
    83   real,dimension(0:nymaxn-1) :: cosf
     102!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
    84103
    85104  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
    86105
    87106  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
     107  integer :: nxx, nyy ! max of nest where we are working
    88108  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
    89109
    … …  
    93113  real,parameter :: const=r_air/ga
    94114  real :: tot_cloud_h
    95   integer :: nxm1, nym1
    96 
    97 !  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
     115
     116!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagn.
    98117
    99118! Loop over all nests
    100119!********************
    101120
    102   do l=1,numbnests
    103     nxm1=nxn(l)-1
    104     nym1=nyn(l)-1
    105 
    106 ! Loop over the whole grid
     121  l_loop: do l=1,numbnests
     122    nxx=nxn(l)-1 ! shorthand
     123    nyy=nyn(l)-1 ! shorthand
     124
     125! Loop over the whole grid (partly implicitly
    107126!*************************
    108 
    109     do jy=0,nyn(l)-1
    110       do ix=0,nxn(l)-1
    111         tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    112              psn(ix,jy,1,n,l))
    113       end do
    114     end do
    115 
    116     pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l)
    117     uvzlev(0:nxm1,0:nym1,1)=0.
    118     wzlev(0:nxm1,0:nym1,1)=0.
    119     rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1))
     127! initialise the automatic arrays
     128
     129  tvold(0:nxx,0:nyy)=tt2n(0:nxx,0:nyy,1,n,l) * &
     130    (1.+0.378*ew( td2n(0:nxx,0:nyy,1,n,l) ) / psn(0:nxx,0:nyy,1,n,l))
     131  pold(0:nxx,0:nyy)=psn(0:nxx,0:nyy,1,n,l)
     132  uvzlev(0:nxx,0:nyy,1)=0.
     133  wzlev(0:nxx,0:nyy,1)=0.
     134  rhohn(0:nxx,0:nyy,1)=pold(0:nxx,0:nyy)/(r_air*tvold(0:nxx,0:nyy))
     135
    120136
    121137! Compute heights of eta levels
    … …  
    123139
    124140    do kz=2,nuvz
    125       pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l)
    126       tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l))
    127       rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1))
    128 
    129       where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2)
    130         uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
    131              &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* &
    132              (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/&
    133              &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1))
     141      pint(:,:)=akz(kz)+bkz(kz)*psn(0:nxx,0:nyy,1,n,l)
     142      tv(0:nxx,0:nyy)=tthn(0:nxx,0:nyy,kz,n,l)* &
     143        (1.+0.608*qvhn(0:nxx,0:nyy,kz,n,l))
     144      rhohn(0:nxx,0:nyy,kz)=pint(0:nxx,0:nyy)/(r_air*tv(0:nxx,0:nyy))
     145
     146      where (abs(tv(0:nxx,0:nyy)-tvold(0:nxx,0:nyy)).gt.0.2)
     147        uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
     148          const*log( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) * &
     149          ( tv(0:nxx,0:nyy) - tvold(0:nxx,0:nyy) ) / &
     150          log( tv(0:nxx,0:nyy)/tvold(0:nxx,0:nyy) )
    134151      elsewhere
    135         uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
    136              &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1)
     152        uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
     153           const*log( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) *  &
     154           tv(0:nxx,0:nyy)
    137155      endwhere
    138156
    139       tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1)
    140       pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1)
     157      tvold(0:nxx,0:nyy)=tv(0:nxx,0:nyy)
     158      pold(0:nxx,0:nyy)=pint(0:nxx,0:nyy)
    141159
    142160    end do
    143161
    144162    do kz=2,nwz-1
    145       wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2.
     163      wzlev(0:nxx,0:nyy,kz)=(uvzlev(0:nxx,0:nyy,kz+1)+uvzlev(0:nxx,0:nyy,kz))/2.
    146164    end do
    147     wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ &
    148          uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1)
    149 
    150 
    151     pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ &
    152          ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- &
    153          (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l)))
     165    wzlev(0:nxx,0:nyy,nwz)=wzlev(0:nxx,0:nyy,nwz-1)+ &
     166         uvzlev(0:nxx,0:nyy,nuvz)-uvzlev(0:nxx,0:nyy,nuvz-1)
     167
     168
     169    pinmconv(0:nxx,0:nyy,1)=(uvzlev(0:nxx,0:nyy,2))/ &
     170      ((aknew(2)+bknew(2)*psn(0:nxx,0:nyy,1,n,l))- &
     171       (aknew(1)+bknew(1)*psn(0:nxx,0:nyy,1,n,l)))
    154172    do kz=2,nz-1
    155       pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ &
    156            ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- &
    157            (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l)))
     173      pinmconv(0:nxx,0:nyy,kz)= &
     174        (uvzlev(0:nxx,0:nyy,kz+1)-uvzlev(0:nxx,0:nyy,kz-1))/ &
     175        ((aknew(kz+1)+bknew(kz+1)*psn(0:nxx,0:nyy,1,n,l))- &
     176         (aknew(kz-1)+bknew(kz-1)*psn(0:nxx,0:nyy,1,n,l)))
    158177    end do
    159     pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ &
    160          ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- &
    161          (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l)))
    162 
    163 ! Levels, where u,v,t and q are given
     178    pinmconv(0:nxx,0:nyy,nz)= &
     179      (uvzlev(0:nxx,0:nyy,nz)-uvzlev(0:nxx,0:nyy,nz-1))/ &
     180      ((aknew(nz)+bknew(nz)*psn(0:nxx,0:nyy,1,n,l))- &
     181       (aknew(nz-1)+bknew(nz-1)*psn(0:nxx,0:nyy,1,n,l)))
     182
     183! Levels where u,v,t and q are given
    164184!************************************
    165185
    166     uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l)
    167     vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l)
    168     ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l)
    169     qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l)
     186    uun(0:nxx,0:nyy,1,n,l)=uuhn(0:nxx,0:nyy,1,l)
     187    vvn(0:nxx,0:nyy,1,n,l)=vvhn(0:nxx,0:nyy,1,l)
     188    ttn(0:nxx,0:nyy,1,n,l)=tthn(0:nxx,0:nyy,1,n,l)
     189    qvn(0:nxx,0:nyy,1,n,l)=qvhn(0:nxx,0:nyy,1,n,l)
    170190    if (readclouds_nest(l)) then
    171       clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l)
    172       if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l)
     191      clwcn(0:nxx,0:nyy,1,n,l)=clwchn(0:nxx,0:nyy,1,n,l)
     192      if (.not. sumclouds_nest(l)) &
     193        ciwcn(0:nxx,0:nyy,1,n,l)=ciwchn(0:nxx,0:nyy,1,n,l)
    173194    end if
    174     pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l)
    175     rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1)
    176 
    177     uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l)
    178     vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l)
    179     ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l)
    180     qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l)
     195    pvn(0:nxx,0:nyy,1,n,l)=pvhn(0:nxx,0:nyy,1,l)
     196    rhon(0:nxx,0:nyy,1,n,l)=rhohn(0:nxx,0:nyy,1)
     197
     198    uun(0:nxx,0:nyy,nz,n,l)=uuhn(0:nxx,0:nyy,nuvz,l)
     199    vvn(0:nxx,0:nyy,nz,n,l)=vvhn(0:nxx,0:nyy,nuvz,l)
     200    ttn(0:nxx,0:nyy,nz,n,l)=tthn(0:nxx,0:nyy,nuvz,n,l)
     201    qvn(0:nxx,0:nyy,nz,n,l)=qvhn(0:nxx,0:nyy,nuvz,n,l)
    181202    if (readclouds_nest(l)) then
    182       clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l)
    183       if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l)
     203      clwcn(0:nxx,0:nyy,nz,n,l)=clwchn(0:nxx,0:nyy,nuvz,n,l)
     204      if (.not. sumclouds_nest(l)) &
     205        ciwcn(0:nxx,0:nyy,nz,n,l)=ciwchn(0:nxx,0:nyy,nuvz,n,l)
    184206    end if
    185     pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l)
    186     rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz)
    187 
     207    pvn(0:nxx,0:nyy,nz,n,l)=pvhn(0:nxx,0:nyy,nuvz,l)
     208    rhon(0:nxx,0:nyy,nz,n,l)=rhohn(0:nxx,0:nyy,nuvz)
    188209
    189210    kmin=2
    190211    idx=kmin
    191     do iz=2,nz-1
    192       do jy=0,nyn(l)-1
    193         do ix=0,nxn(l)-1
    194           if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
     212    iz_loop: do iz=2,nz-1
     213
     214      do jy=0,nyy
     215        do ix=0,nxx
     216          if( height(iz).gt.uvzlev(ix,jy,nuvz)) then
     217
    195218            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    196219            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    … …  
    200223            if (readclouds_nest(l)) then
    201224              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
    202               if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
     225              if (.not. sumclouds_nest(l))  &
     226                  ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
    203227            end if
    204228!hg
    205229            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
    206230            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
     231
    207232          else
    208             innuvz: do kz=idx(ix,jy),nuvz
    209               if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    210                    (height(iz).le.uvzlev(ix,jy,kz))) then
     233
     234            kz_loop: do kz=idx(ix,jy),nuvz
     235              if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
     236                  .and. height(iz).le.uvzlev(ix,jy,kz)) then
    211237                idx(ix,jy)=kz
    212                 exit innuvz
     238                exit kz_loop
    213239              endif
    214             enddo innuvz
     240            enddo kz_loop
     241           
    215242          endif
    216243        enddo
    217244      enddo
    218       do jy=0,nyn(l)-1
    219         do ix=0,nxn(l)-1
     245
     246      do jy=0,nyy
     247        do ix=0,nxx
    220248          if(height(iz).le.uvzlev(ix,jy,nuvz)) then
    221249            kz=idx(ix,jy)
    … …  
    231259!hg adding the cloud water
    232260            if (readclouds_nest(l)) then
    233               clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz
    234               if (.not.sumclouds_nest(l)) &
    235                    &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
     261              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2 + &
     262                clwchn(ix,jy,kz,n,l)*dz1)/dz
     263              if (.not. sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)= &
     264                  (ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
    236265            end if
    237266!hg
    238267            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
    239268            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
     269           
    240270          endif
    241271        enddo
    242272      enddo
    243     enddo
     273
     274    enddo iz_loop
    244275
    245276
    … …  
    282313!*************************
    283314
    284     wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1)
    285     wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz)
     315    wwn(0:nxx,0:nyy,1,n,l)=wwhn(0:nxx,0:nyy,1,l)*pinmconv(0:nxx,0:nyy,1)
     316    wwn(0:nxx,0:nyy,nz,n,l)=wwhn(0:nxx,0:nyy,nwz,l)*pinmconv(0:nxx,0:nyy,nz)
    286317    kmin=2
    287318    idx=kmin
    288     do iz=2,nz
    289       do jy=0,nym1
    290         do ix=0,nxm1
    291           inn:         do kz=idx(ix,jy),nwz
    292             if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
    293                 height(iz).le.wzlev(ix,jy,kz)) then
     319    iz_loop2: do iz=2,nz
     320      do jy=0,nyy
     321        do ix=0,nxx
     322          kz_loop2:   do kz=idx(ix,jy),nwz
     323            if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
     324                .and. height(iz).le.wzlev(ix,jy,kz)) then
    294325              idx(ix,jy)=kz
    295               exit inn
     326              exit kz_loop2
    296327            endif
    297           enddo inn
     328          enddo kz_loop2
    298329        enddo
    299330      enddo
    300       do jy=0,nym1
    301         do ix=0,nxm1
     331      do jy=0,nyy
     332        do ix=0,nxx
    302333          kz=idx(ix,jy)
    303334          dz1=height(iz)-wzlev(ix,jy,kz-1)
    … …  
    308339        enddo
    309340      enddo
    310     enddo
     341    enddo iz_loop2
    311342
    312343!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
    … …  
    332363!*************************************************
    333364
    334     drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ &
    335          (height(2)-height(1))
     365    drhodzn(0:nxx,0:nyy,1,n,l)= &
     366      ( rhon(0:nxx,0:nyy,2,n,l)-rhon(0:nxx,0:nyy,1,n,l) )/(height(2)-height(1))
    336367    do kz=2,nz-1
    337       drhodzn(0:nxm1,0:nym1,kz,n,l)=(rhon(0:nxm1,0:nym1,kz+1,n,l)-rhon(0:nxm1,0:nym1,kz-1,n,l))/ &
    338            (height(kz+1)-height(kz-1))
     368      drhodzn(0:nxx,0:nyy,kz,n,l)= &
     369        ( rhon(0:nxx,0:nyy,kz+1,n,l)-rhon(0:nxx,0:nyy,kz-1,n,l) ) / &
     370        ( height(kz+1)-height(kz-1) )
    339371    end do
    340     drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l)
     372    drhodzn(0:nxx,0:nyy,nz,n,l)=drhodzn(0:nxx,0:nyy,nz-1,n,l)
    341373
    342374
    … …  
    360392!****************************************************************
    361393
    362     do jy=1,nyn(l)-2
     394    do jy=1,nyy-1
    363395!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    364       cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     396      cosf(jy) = 1. / cos( ( real(jy)*dyn(l) + ylat0n(l) ) * pi180 )
    365397    end do
    366398
    367399    kmin=2
    368400    idx=kmin
    369     do iz=2,nz-1
    370       do jy=1,nyn(l)-2
    371         do ix=1,nxn(l)-2
    372 
    373           inneta: do kz=idx(ix,jy),nz
    374             if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    375                 (height(iz).le.uvzlev(ix,jy,kz))) then
     401    iz_loop3: do iz=2,nz-1
     402      do jy=1,nyy-1
     403        do ix=1,nxx-1
     404
     405          kz_loop3: do kz=idx(ix,jy),nz
     406            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)) &
     407                .and. (height(iz).le.uvzlev(ix,jy,kz))) then
    376408              idx(ix,jy)=kz
    377               exit inneta
     409              exit kz_loop3
    378410            endif
    379           enddo inneta
     411          enddo kz_loop3
    380412        enddo
    381413      enddo
    382414
    383       do jy=1,nyn(l)-2
    384         do ix=1,nxn(l)-2
     415      do jy=1,nyy-1
     416        do ix=1,nxx-1
    385417          kz=idx(ix,jy)
    386418          dz1=height(iz)-uvzlev(ix,jy,kz-1)
    … …  
    400432          dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    401433
    402           wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+&
    403                &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l))
     434          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l) + &
     435            ( dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy) + &
     436              dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l) )
    404437
    405438        end do
    406439
    407440      end do
    408     end do
     441    end do iz_loop3
    409442
    410443
    … …  
    455488
    456489
    457 !*********************************************************************************** 
    458     if (readclouds_nest(l)) then !HG METHOD
    459 ! The method is loops all grids vertically and constructs the 3D matrix for clouds
    460 ! Cloud top and cloud bottom gid cells are assigned as well as the total column
    461 ! cloud water. For precipitating grids, the type and whether it is in or below
     490!******************************************************************************
     491    if (readclouds_nest(l)) then !HG METHOD
     492   
     493! Loops over all grid cells vertically and construct the 3D matrix for clouds
     494! Cloud top and cloud bottom grid cells are assigned as well as the total column
     495! cloud water. For precipitating columns, the type and whether it is in or below
    462496! cloud scavenging are assigned with numbers 2-5 (following the old metod).
    463 ! Distinction is done for lsp and convp though they are treated the same in regards
    464 ! to scavenging. Also clouds that are not precipitating are defined which may be
    465 ! to include future cloud processing by non-precipitating-clouds.
    466 !***********************************************************************************
    467       write(*,*) 'Nested ECMWF fields: using cloud water'
    468       clwn(0:nxm1,0:nym1,:,n,l)=0.0
    469 !    icloud_stats(0:nxm1,0:nym1,:,n)=0.0
    470       ctwcn(0:nxm1,0:nym1,n,l)=0.0
    471       cloudsn(0:nxm1,0:nym1,:,n,l)=0
     497! A distinction is made between lsp and convp though they are treated the equally
     498! with regard to scavenging. Also, clouds that are not precipitating are defined which
     499! may be used in the future for cloud processing by non-precipitating-clouds.
     500!*******************************************************************************
     501
     502!PS      write(*,*) 'Nested ECMWF fields: using cloud water'
     503      clwn(0:nxx,0:nyy,:,n,l)=0.
     504!    icloud_stats(0:nxx,0:nyy,:,n)=0.
     505      ctwcn(0:nxx,0:nyy,n,l)=0.
     506      cloudsn(0:nxx,0:nyy,:,n,l)=0
    472507! If water/ice are read separately into clwc and ciwc, store sum in clwcn
    473       if (.not.sumclouds_nest(l)) then
    474         clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l)
     508      if (.not. sumclouds_nest(l)) then
     509        clwcn(0:nxx,0:nyy,:,n,l) = clwcn(0:nxx,0:nyy,:,n,l) + &
     510                                   ciwcn(0:nxx,0:nyy,:,n,l)
    475511      end if
    476       do jy=0,nym1
    477         do ix=0,nxm1
     512      do jy=0,nyy
     513        do ix=0,nxx
    478514          lsp=lsprecn(ix,jy,1,n,l)
    479515          convp=convprecn(ix,jy,1,n,l)
    … …  
    481517          tot_cloud_h=0
    482518! Find clouds in the vertical
     519!! Note PS: bad loop order.
    483520          do kz=1, nz-1 !go from top to bottom
    484521            if (clwcn(ix,jy,kz,n,l).gt.0) then     
    485 ! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
    486               clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz))
     522! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
     523              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l)) * &
     524                (height(kz+1)-height(kz))
    487525              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
    488526              ctwcn(ix,jy,n,l) = ctwcn(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
    … …  
    499537
    500538! If Precipitation. Define removal type in the vertical
    501           if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     539        if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
     540!! Note PS: such hardcoded limits would better be parameters defined in par_mod
    502541
    503542            do kz=nz,1,-1 !go Bottom up!
    504               if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud
    505                 cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
    506                 cloudsn(ix,jy,kz,n,l)=1                               ! is a cloud
     543!! Note PS: bad loop order
     544             if (clwn(ix,jy,kz,n,l) .gt. 0.) then ! is in cloud
     545               cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
     546               cloudsn(ix,jy,kz,n,l)=1                             ! is a cloud
    507547                if (lsp.ge.convp) then
    508                   cloudsn(ix,jy,kz,n,l)=3                            ! lsp in-cloud
     548                  cloudsn(ix,jy,kz,n,l)=3                        ! lsp in-cloud
    509549                else
    510                   cloudsn(ix,jy,kz,n,l)=2                             ! convp in-cloud
    511                 endif                                              ! convective or large scale
    512               else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
    513                 if (lsp.ge.convp) then
    514                   cloudsn(ix,jy,kz,n,l)=5                             ! lsp dominated washout
     550                  cloudsn(ix,jy,kz,n,l)=2                      ! convp in-cloud
     551                endif                               ! convective or large scale
     552              elseif ( clwn(ix,jy,kz,n,l) .le. 0. .and. &
     553                  cloudh_min .ge. height(kz) ) then ! is below cloud
     554                if (lsp .ge. convp) then
     555                  cloudsn(ix,jy,kz,n,l)=5              ! lsp dominated washout
    515556                else
    516                   cloudsn(ix,jy,kz,n,l)=4                             ! convp dominated washout
    517                 endif                                              ! convective or large scale
     557                  cloudsn(ix,jy,kz,n,l)=4            ! convp dominated washout
     558                endif                             ! convective or large scale
    518559              endif
    519560
    520               if (height(kz).ge. 19000) then                        ! set a max height for removal
     561              if (height(kz).ge. 19000) then ! set a max height for removal
    521562                cloudsn(ix,jy,kz,n,l)=0
    522               endif !clw>0
    523             end do !nz
     563              endif ! clw>0
     564            end do ! kz
    524565          endif ! precipitation
    525566        end do
    526567      end do
     568     
    527569!********************************************************************
    528570    else ! old method:
    529571!********************************************************************
    530       write(*,*) 'Nested fields: using cloud water from Parameterization'
    531       do jy=0,nyn(l)-1
    532         do ix=0,nxn(l)-1
     572 
     573!PS      write(*,*) 'Nested fields: using cloud water from Parameterization'
     574      do jy=0,nyy
     575        do ix=0,nxx
    533576          rain_cloud_above(ix,jy)=0
    534577          lsp=lsprecn(ix,jy,1,n,l)
    … …  
    536579          cloudshn(ix,jy,n,l)=0
    537580          do kz_inv=1,nz-1
     581!! Note PS: bad loop order.
    538582            kz=nz-kz_inv+1
    539583            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
    540584            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
    541585            cloudsn(ix,jy,kz,n,l)=0
    542             if (rh.gt.0.8) then ! in cloud
    543               if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
     586
     587            if (rh .gt. 0.8) then ! in cloud
     588!! Note PS: such hardcoded limits would better be parameters defined in par_mod
     589
     590              if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
     591!! Note PS: such hardcoded limits would better be parameters defined in par_mod
    544592                rain_cloud_above(ix,jy)=1
    545593                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
    … …  
    553601                cloudsn(ix,jy,kz,n,l)=1 ! cloud
    554602              endif
     603
    555604            else ! no cloud
     605
    556606              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
    557607                if (lsp.ge.convp) then
    … …  
    561611                endif
    562612              endif
     613
    563614            endif
    564           end do
    565         end do
    566       end do
    567     end if
    568 
    569   end do ! end loop over nests
     615          end do ! kz
     616         
     617        end do ! ix
     618      end do ! jy
     619     
     620    end if! old method:
     621
     622  end do l_loop! end loop over nests
    570623
    571624end subroutine verttransform_nests
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG