Changeset db712a8 in flexpart.git for src/verttransform_nests.f90


Ignore:
Timestamp:
Mar 3, 2016, 12:34:56 PM (8 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
38b7917
Parents:
b0434e1
Message:

Completed handling of nested wind fields with cloud water (for wet deposition).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/verttransform_nests.f90

    rb0434e1 rdb712a8  
    2121
    2222subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
    23   !                            i   i    i    i   i
    24   !*****************************************************************************
    25   !                                                                            *
    26   !     This subroutine transforms temperature, dew point temperature and      *
    27   !     wind components from eta to meter coordinates.                         *
    28   !     The vertical wind component is transformed from Pa/s to m/s using      *
    29   !     the conversion factor pinmconv.                                        *
    30   !     In addition, this routine calculates vertical density gradients        *
    31   !     needed for the parameterization of the turbulent velocities.           *
    32   !     It is similar to verttransform, but makes the transformations for      *
    33   !     the nested grids.                                                      *
    34   !                                                                            *
    35   !     Author: A. Stohl, G. Wotawa                                            *
    36   !                                                                            *
    37   !     12 August 1996                                                         *
    38   !     Update: 16 January 1998                                                *
    39   !                                                                            *
    40   !     Major update: 17 February 1999                                         *
    41   !     by G. Wotawa                                                           *
    42   !                                                                            *
    43   !     - Vertical levels for u, v and w are put together                      *
    44   !     - Slope correction for vertical velocity: Modification of calculation  *
    45   !       procedure                                                            *
    46   !                                                                            *
    47   !*****************************************************************************
    48   !  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
    49   !   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   !*****************************************************************************
    54   !                                                                            *
    55   ! Variables:                                                                 *
    56   ! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
    57   ! uun                             wind components in x-direction [m/s]       *
    58   ! vvn                             wind components in y-direction [m/s]       *
    59   ! wwn                             wind components in z-direction [deltaeta/s]*
    60   ! ttn                             temperature [K]                            *
    61   ! pvn                             potential vorticity (pvu)                  *
    62   ! psn                             surface pressure [Pa]                      *
    63   !                                                                            *
    64   !*****************************************************************************
     23!                            i   i    i    i   i
     24!*****************************************************************************
     25!                                                                            *
     26!     This subroutine transforms temperature, dew point temperature and      *
     27!     wind components from eta to meter coordinates.                         *
     28!     The vertical wind component is transformed from Pa/s to m/s using      *
     29!     the conversion factor pinmconv.                                        *
     30!     In addition, this routine calculates vertical density gradients        *
     31!     needed for the parameterization of the turbulent velocities.           *
     32!     It is similar to verttransform, but makes the transformations for      *
     33!     the nested grids.                                                      *
     34!                                                                            *
     35!     Author: A. Stohl, G. Wotawa                                            *
     36!                                                                            *
     37!     12 August 1996                                                         *
     38!     Update: 16 January 1998                                                *
     39!                                                                            *
     40!     Major update: 17 February 1999                                         *
     41!     by G. Wotawa                                                           *
     42!                                                                            *
     43!     - Vertical levels for u, v and w are put together                      *
     44!     - Slope correction for vertical velocity: Modification of calculation  *
     45!       procedure                                                            *
     46!                                                                            *
     47!*****************************************************************************
     48!  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
     49!   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!*****************************************************************************
     54! ESO, 2016
     55! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than
     56!  the actual field dimensions
     57!*****************************************************************************
     58!                                                                            *
     59! Variables:                                                                 *
     60! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
     61! uun                             wind components in x-direction [m/s]       *
     62! vvn                             wind components in y-direction [m/s]       *
     63! wwn                             wind components in z-direction [deltaeta/s]*
     64! ttn                             temperature [K]                            *
     65! pvn                             potential vorticity (pvu)                  *
     66! psn                             surface pressure [Pa]                      *
     67!                                                                            *
     68!*****************************************************************************
    6569
    6670  use par_mod
    … …  
    6973  implicit none
    7074
    71   integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp
    72   integer :: rain_cloud_above,kz_inv
    73   real :: f_qvsat,pressure,rh,lsp,convp
    74   real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
    75   real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
    76   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
     75  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn
     76  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
     77
     78  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
     79  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
     80  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
     81  real,dimension(0:nymaxn-1) :: cosf
     82
     83  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
     84
     85  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
     86  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
     87
     88  real :: ew,dz1,dz2,dz
    7789  real :: dzdx,dzdy
    7890  real :: dzdx1,dzdx2,dzdy1,dzdy2
    79   real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    80   real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    81   real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    82   real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
    8391  real,parameter :: const=r_air/ga
    84 
    85 
    86   ! Loop over all nests
    87   !********************
     92  real :: tot_cloud_h
     93  integer :: nxm1, nym1
     94
     95!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
     96
     97! Loop over all nests
     98!********************
    8899
    89100  do l=1,numbnests
    90 
    91   ! Loop over the whole grid
    92   !*************************
    93 
    94   do jy=0,nyn(l)-1
    95     do ix=0,nxn(l)-1
    96 
    97       tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    98       psn(ix,jy,1,n,l))
    99       pold=psn(ix,jy,1,n,l)
    100       uvwzlev(ix,jy,1)=0.
    101       wzlev(1)=0.
    102       rhoh(1)=pold/(r_air*tvold)
    103 
    104 
    105   ! Compute heights of eta levels
    106   !******************************
    107 
    108       do kz=2,nuvz
    109         pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)
    110         tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
    111         rhoh(kz)=pint/(r_air*tv)
    112 
    113         if (abs(tv-tvold).gt.0.2) then
    114           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    115           (tv-tvold)/log(tv/tvold)
    116         else
    117           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    118         endif
    119 
    120         tvold=tv
    121         pold=pint
     101    nxm1=nxn(l)-1
     102    nym1=nyn(l)-1
     103
     104! Loop over the whole grid
     105!*************************
     106
     107    do jy=0,nyn(l)-1
     108      do ix=0,nxn(l)-1
     109        tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew*(td2n(ix,jy,1,n,l))/ &
     110             psn(ix,jy,1,n,l))
    122111      end do
    123 
    124 
    125       do kz=2,nwz-1
    126         wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
    127       end do
    128       wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    129 
    130   ! pinmconv=(h2-h1)/(p2-p1)
    131 
    132       pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    133       ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    134       (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    135       do kz=2,nz-1
    136         pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    137         ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
    138         (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
    139       end do
    140       pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    141       ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    142       (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    143 
    144 
    145   ! Levels, where u,v,t and q are given
    146   !************************************
    147 
    148       uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
    149       vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
    150       ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
    151       qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
    152       pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
    153       rhon(ix,jy,1,n,l)=rhoh(1)
    154       uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
    155       vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
    156       ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
    157       qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
    158       pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
    159       rhon(ix,jy,nz,n,l)=rhoh(nuvz)
    160       kmin=2
    161       do iz=2,nz-1
    162         do kz=kmin,nuvz
    163           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     112    end do
     113
     114    pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l)
     115    uvzlev(0:nxm1,0:nym1,1)=0.
     116    wzlev(0:nxm1,0:nym1,1)=0.
     117    rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1))
     118
     119! Compute heights of eta levels
     120!******************************
     121
     122    do kz=2,nuvz
     123      pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l)
     124      tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l))
     125      rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1))
     126
     127      where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2)
     128        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
     129             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* &
     130             (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/&
     131             &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1))
     132      elsewhere
     133        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
     134             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1)
     135      endwhere
     136
     137      tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1)
     138      pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1)
     139
     140    end do
     141
     142    do kz=2,nwz-1
     143      wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2.
     144    end do
     145    wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ &
     146         uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1)
     147
     148
     149    pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ &
     150         ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- &
     151         (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l)))
     152    do kz=2,nz-1
     153      pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ &
     154           ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- &
     155           (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l)))
     156    end do
     157    pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ &
     158         ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- &
     159         (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l)))
     160
     161! Levels, where u,v,t and q are given
     162!************************************
     163
     164    uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l)
     165    vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l)
     166    ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l)
     167    qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l)
     168    if (readclouds_nest(l)) then
     169      clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l)
     170      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l)
     171    end if
     172    pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l)
     173    rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1)
     174
     175    uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l)
     176    vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l)
     177    ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l)
     178    qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l)
     179    if (readclouds_nest(l)) then
     180      clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l)
     181      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l)
     182    end if
     183    pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l)
     184    rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz)
     185
     186
     187    kmin=2
     188    idx=kmin
     189    do iz=2,nz-1
     190      do jy=0,nyn(l)-1
     191        do ix=0,nxn(l)-1
     192          if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
    164193            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    165194            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    166195            ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
    167196            qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
     197!hg adding the cloud water
     198            if (readclouds_nest(l)) then
     199              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
     200              if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
     201            end if
     202!hg
    168203            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
    169204            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
    170             goto 30
     205          else
     206            innuvz: do kz=idx(ix,jy),nuvz
     207              if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
     208                   (height(iz).le.uvzlev(ix,jy,kz))) then
     209                idx(ix,jy)=kz
     210                exit innuvz
     211              endif
     212            enddo innuvz
    171213          endif
    172           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    173           (height(iz).le.uvwzlev(ix,jy,kz))) then
    174            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    175            dz2=uvwzlev(ix,jy,kz)-height(iz)
    176            dz=dz1+dz2
    177            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    178            uuhn(ix,jy,kz,l)*dz1)/dz
    179            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    180            vvhn(ix,jy,kz,l)*dz1)/dz
    181            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    182            tthn(ix,jy,kz,n,l)*dz1)/dz
    183            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    184            qvhn(ix,jy,kz,n,l)*dz1)/dz
    185            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    186            pvhn(ix,jy,kz,l)*dz1)/dz
    187            rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    188            kmin=kz
    189            goto 30
     214        enddo
     215      enddo
     216      do jy=0,nyn(l)-1
     217        do ix=0,nxn(l)-1
     218          if(height(iz).le.uvzlev(ix,jy,nuvz)) then
     219            kz=idx(ix,jy)
     220            dz1=height(iz)-uvzlev(ix,jy,kz-1)
     221            dz2=uvzlev(ix,jy,kz)-height(iz)
     222            dz=dz1+dz2
     223            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz
     224            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz
     225            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
     226                 +tthn(ix,jy,kz,n,l)*dz1)/dz
     227            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
     228                 +qvhn(ix,jy,kz,n,l)*dz1)/dz
     229!hg adding the cloud water
     230            if (readclouds_nest(l)) then
     231              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz
     232              if (.not.sumclouds_nest(l)) &
     233                   &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
     234            end if
     235!hg
     236            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
     237            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
    190238          endif
     239        enddo
     240      enddo
     241    enddo
     242
     243
     244
     245!         do kz=kmin,nuvz
     246!           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     247!             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
     248!             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
     249!             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
     250!             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
     251!             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
     252!             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
     253!             goto 30
     254!           endif
     255!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     256!           (height(iz).le.uvwzlev(ix,jy,kz))) then
     257!            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     258!            dz2=uvwzlev(ix,jy,kz)-height(iz)
     259!            dz=dz1+dz2
     260!            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
     261!            uuhn(ix,jy,kz,l)*dz1)/dz
     262!            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
     263!            vvhn(ix,jy,kz,l)*dz1)/dz
     264!            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
     265!            tthn(ix,jy,kz,n,l)*dz1)/dz
     266!            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
     267!            qvhn(ix,jy,kz,n,l)*dz1)/dz
     268!            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
     269!            pvhn(ix,jy,kz,l)*dz1)/dz
     270!            rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz
     271!            kmin=kz
     272!            goto 30
     273!           endif
     274!         end do
     275! 30      continue
     276!       end do
     277
     278
     279! Levels, where w is given
     280!*************************
     281
     282    wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1)
     283    wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz)
     284    kmin=2
     285    idx=kmin
     286    do iz=2,nz
     287      do jy=0,nym1
     288        do ix=0,nxm1
     289          inn:         do kz=idx(ix,jy),nwz
     290            if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
     291                 height(iz).le.wzlev(ix,jy,kz)) then
     292              idx(ix,jy)=kz
     293              exit inn
     294            endif
     295          enddo inn
     296        enddo
     297      enddo
     298      do jy=0,nym1
     299        do ix=0,nxm1
     300          kz=idx(ix,jy)
     301          dz1=height(iz)-wzlev(ix,jy,kz-1)
     302          dz2=wzlev(ix,jy,kz)-height(iz)
     303          dz=dz1+dz2
     304          wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 &
     305               +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz
     306        enddo
     307      enddo
     308    enddo
     309
     310!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
     311!       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
     312!       kmin=2
     313!       do iz=2,nz
     314!         do kz=kmin,nwz
     315!           if ((height(iz).gt.wzlev(kz-1)).and. &
     316!           (height(iz).le.wzlev(kz))) then
     317!             dz1=height(iz)-wzlev(kz-1)
     318!             dz2=wzlev(kz)-height(iz)
     319!             dz=dz1+dz2
     320!             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
     321!             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
     322!             kmin=kz
     323!             goto 40
     324!           endif
     325!         end do
     326! 40      continue
     327!       end do
     328
     329! Compute density gradients at intermediate levels
     330!*************************************************
     331
     332    drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ &
     333         (height(2)-height(1))
     334    do kz=2,nz-1
     335      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))/ &
     336           (height(kz+1)-height(kz-1))
     337    end do
     338    drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l)
     339
     340
     341! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
     342! (height(2)-height(1))
     343! do kz=2,nz-1
     344!   drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
     345!   rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     346! end do
     347! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
     348
     349
     350
     351!   end do
     352! end do
     353
     354
     355!****************************************************************
     356! Compute slope of eta levels in windward direction and resulting
     357! vertical wind correction
     358!****************************************************************
     359
     360    do jy=1,nyn(l)-2
     361!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     362      cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     363    end do
     364
     365    kmin=2
     366    idx=kmin
     367    do iz=2,nz-1
     368      do jy=1,nyn(l)-2
     369        do ix=1,nxn(l)-2
     370
     371          inneta: do kz=idx(ix,jy),nz
     372            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
     373                 (height(iz).le.uvzlev(ix,jy,kz))) then
     374              idx(ix,jy)=kz
     375              exit inneta
     376            endif
     377          enddo inneta
     378        enddo
     379      enddo
     380
     381      do jy=1,nyn(l)-2
     382        do ix=1,nxn(l)-2
     383          kz=idx(ix,jy)
     384          dz1=height(iz)-uvzlev(ix,jy,kz-1)
     385          dz2=uvzlev(ix,jy,kz)-height(iz)
     386          dz=dz1+dz2
     387          ix1=ix-1
     388          jy1=jy-1
     389          ixp=ix+1
     390          jyp=jy+1
     391
     392          dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
     393          dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
     394          dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
     395
     396          dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
     397          dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
     398          dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
     399
     400          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+&
     401               &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l))
     402
    191403        end do
    192 30      continue
     404
    193405      end do
    194 
    195 
    196   ! Levels, where w is given
    197   !*************************
    198 
    199       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
    200       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
    201       kmin=2
    202       do iz=2,nz
    203         do kz=kmin,nwz
    204           if ((height(iz).gt.wzlev(kz-1)).and. &
    205           (height(iz).le.wzlev(kz))) then
    206             dz1=height(iz)-wzlev(kz-1)
    207             dz2=wzlev(kz)-height(iz)
    208             dz=dz1+dz2
    209             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
    210             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
    211             kmin=kz
    212             goto 40
    213           endif
     406    end do
     407
     408
     409!   do jy=1,nyn(l)-2
     410!     do ix=1,nxn(l)-2
     411!       kmin=2
     412!       do iz=2,nz-1
     413
     414!         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy)
     415!         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
     416
     417!         do kz=kmin,nz
     418!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     419!                (height(iz).le.uvwzlev(ix,jy,kz))) then
     420!             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     421!             dz2=uvwzlev(ix,jy,kz)-height(iz)
     422!             dz=dz1+dz2
     423!             kl=kz-1
     424!             klp=kz
     425!             kmin=kz
     426!             goto 47
     427!           endif
     428!         end do
     429
     430! 47      ix1=ix-1
     431!         jy1=jy-1
     432!         ixp=ix+1
     433!         jyp=jy+1
     434
     435!         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
     436!         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
     437!         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
     438
     439!         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
     440!         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
     441!         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
     442
     443!         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
     444
     445!       end do
     446
     447!     end do
     448!   end do
     449
     450
     451!write (*,*) 'initializing nested cloudsn, n:',n
     452!   create a cloud and rainout/washout field, cloudsn occur where rh>80%
     453
     454
     455!*********************************************************************************** 
     456    if (readclouds_nest(l)) then !HG METHOD
     457! The method is loops all grids vertically and constructs the 3D matrix for clouds
     458! Cloud top and cloud bottom gid cells are assigned as well as the total column
     459! cloud water. For precipitating grids, the type and whether it is in or below
     460! cloud scavenging are assigned with numbers 2-5 (following the old metod).
     461! Distinction is done for lsp and convp though they are treated the same in regards
     462! to scavenging. Also clouds that are not precipitating are defined which may be
     463! to include future cloud processing by non-precipitating-clouds.
     464!***********************************************************************************
     465      write(*,*) 'Nested ECMWF fields: using cloud water'
     466      clwn(0:nxm1,0:nym1,:,n,l)=0.0
     467!    icloud_stats(0:nxm1,0:nym1,:,n)=0.0
     468      clw4n(0:nxm1,0:nym1,n,l)=0.0
     469      cloudsn(0:nxm1,0:nym1,:,n,l)=0
     470! If water/ice are read separately into clwc and ciwc, store sum in clwcn
     471      if (.not.sumclouds_nest(l)) then
     472        clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l)
     473      end if
     474      do jy=0,nym1
     475        do ix=0,nxm1
     476          lsp=lsprecn(ix,jy,1,n,l)
     477          convp=convprecn(ix,jy,1,n,l)
     478          prec=lsp+convp
     479          tot_cloud_h=0
     480! Find clouds in the vertical
     481          do kz=1, nz-1 !go from top to bottom
     482            if (clwcn(ix,jy,kz,n,l).gt.0) then     
     483! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
     484              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz))
     485              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
     486              clw4n(ix,jy,n,l) = clw4n(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
     487!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
     488!           icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
     489              cloudh_min=min(height(kz+1),height(kz))
     490!ZHG 2015 extra for testing
     491!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
     492!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
     493!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
     494!ZHG
     495            endif
     496          end do
     497
     498! If Precipitation. Define removal type in the vertical
     499          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     500
     501            do kz=nz,1,-1 !go Bottom up!
     502              if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud
     503                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
     504                cloudsn(ix,jy,kz,n,l)=1                               ! is a cloud
     505                if (lsp.ge.convp) then
     506                  cloudsn(ix,jy,kz,n,l)=3                            ! lsp in-cloud
     507                else
     508                  cloudsn(ix,jy,kz,n,l)=2                             ! convp in-cloud
     509                endif                                              ! convective or large scale
     510              else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
     511                if (lsp.ge.convp) then
     512                  cloudsn(ix,jy,kz,n,l)=5                             ! lsp dominated washout
     513                else
     514                  cloudsn(ix,jy,kz,n,l)=4                             ! convp dominated washout
     515                endif                                              ! convective or large scale
     516              endif
     517
     518              if (height(kz).ge. 19000) then                        ! set a max height for removal
     519                cloudsn(ix,jy,kz,n,l)=0
     520              endif !clw>0
     521            end do !nz
     522          endif ! precipitation
    214523        end do
    215 40      continue
    216524      end do
    217 
    218   ! Compute density gradients at intermediate levels
    219   !*************************************************
    220 
    221       drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    222       (height(2)-height(1))
    223       do kz=2,nz-1
    224         drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    225         rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     525!********************************************************************
     526    else ! old method:
     527!********************************************************************
     528      write(*,*) 'Nested fields: using cloud water from Parameterization'
     529      do jy=0,nyn(l)-1
     530        do ix=0,nxn(l)-1
     531          rain_cloud_above(ix,jy)=0
     532          lsp=lsprecn(ix,jy,1,n,l)
     533          convp=convprecn(ix,jy,1,n,l)
     534          cloudshn(ix,jy,n,l)=0
     535          do kz_inv=1,nz-1
     536            kz=nz-kz_inv+1
     537            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
     538            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
     539            cloudsn(ix,jy,kz,n,l)=0
     540            if (rh.gt.0.8) then ! in cloud
     541              if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
     542                rain_cloud_above(ix,jy)=1
     543                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
     544                     height(kz)-height(kz-1)
     545                if (lsp.ge.convp) then
     546                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
     547                else
     548                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
     549                endif
     550              else ! no precipitation
     551                cloudsn(ix,jy,kz,n,l)=1 ! cloud
     552              endif
     553            else ! no cloud
     554              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
     555                if (lsp.ge.convp) then
     556                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
     557                else
     558                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
     559                endif
     560              endif
     561            endif
     562          end do
     563        end do
    226564      end do
    227       drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
    228 
    229     end do
    230   end do
    231 
    232 
    233   !****************************************************************
    234   ! Compute slope of eta levels in windward direction and resulting
    235   ! vertical wind correction
    236   !****************************************************************
    237 
    238   do jy=1,nyn(l)-2
    239     cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    240     do ix=1,nxn(l)-2
    241 
    242       kmin=2
    243       do iz=2,nz-1
    244 
    245         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
    246         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    247 
    248         do kz=kmin,nz
    249           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    250                (height(iz).le.uvwzlev(ix,jy,kz))) then
    251             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    252             dz2=uvwzlev(ix,jy,kz)-height(iz)
    253             dz=dz1+dz2
    254             kl=kz-1
    255             klp=kz
    256             kmin=kz
    257             goto 47
    258           endif
    259         end do
    260 
    261 47      ix1=ix-1
    262         jy1=jy-1
    263         ixp=ix+1
    264         jyp=jy+1
    265 
    266         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
    267         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
    268         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
    269 
    270         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
    271         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
    272         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    273 
    274         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
    275 
    276       end do
    277 
    278     end do
    279   end do
    280 
    281 
    282   !write (*,*) 'initializing nested cloudsn, n:',n
    283   !   create a cloud and rainout/washout field, cloudsn occur where rh>80%
    284   write(*,*) 'Nested fields: using cloud water from Parameterization'
    285   do jy=0,nyn(l)-1
    286     do ix=0,nxn(l)-1
    287       rain_cloud_above=0
    288       lsp=lsprecn(ix,jy,1,n,l)
    289       convp=convprecn(ix,jy,1,n,l)
    290       cloudsnh(ix,jy,n,l)=0
    291       do kz_inv=1,nz-1
    292          kz=nz-kz_inv+1
    293          pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
    294          rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
    295          cloudsn(ix,jy,kz,n,l)=0
    296          if (rh.gt.0.8) then ! in cloud
    297             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
    298                rain_cloud_above=1
    299                cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    300                height(kz)-height(kz-1)
    301                if (lsp.ge.convp) then
    302                   cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
    303                else
    304                   cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
    305                endif
    306             else ! no precipitation
    307                   cloudsn(ix,jy,kz,n,l)=1 ! cloud
    308             endif
    309          else ! no cloud
    310             if (rain_cloud_above.eq.1) then ! scavenging
    311                if (lsp.ge.convp) then
    312                   cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
    313                else
    314                   cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
    315                endif
    316             endif
    317          endif
    318       end do
    319     end do
    320   end do
    321 
    322   end do
     565    end if
     566
     567  end do ! end loop over nests
    323568
    324569end subroutine verttransform_nests
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG