Changeset 1a8fbee in flexpart.git for src/verttransform_nests.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_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