Changes in / [b255cd0:ca350ba] in flexpart.git


Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/verttransform.f90

    r4d45639 r0f20c31  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23 !                         i  i   i   i   i
     23!                      i  i   i   i   i
    2424!*****************************************************************************
    2525!                                                                            *
    … …  
    6767  use com_mod
    6868  use cmapf_mod, only: cc2gll
    69 ! eso TODO:
    70 ! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below)
    71 ! as this routine should not be dependent on MPI
    72 !  use mpi_mod
    73 ! :TODO
    7469
    7570  implicit none
    7671
    7772  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    78   integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1)
    79   real :: f_qvsat,pressure,cosf(0:nymax-1)
    80   real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec
    81   real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax)
    82   real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax)
    83   real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1)
    84   real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi
     73  integer :: rain_cloud_above,kz_inv
     74! integer :: icloudtop !PS
     75  real :: f_qvsat,pressure
     76! real :: rh,lsp,convp
     77  real :: rh,lsp,convp,prec,rhmin,precmin
     78  real :: rhoh(nuvzmax),pinmconv(nzmax)
     79  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    8580  real :: xlon,ylat,xlonr,dzdx,dzdy
    86   real :: dzdx1,dzdx2,dzdy1,dzdy2
     81  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    8782  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8883  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    9085  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
    9186  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    92   real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)
     87  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
     88  logical lconvectprec
    9389  real,parameter :: const=r_air/ga
    94 !  integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2
    9590  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9691
    … …  
    108103! If verttransform is called the first time, initialize heights of the   *
    109104! z levels in meter. The heights are the heights of model levels, where  *
    110 ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
    111 ! the vertical resolution in the z system is doubled. As reference point,*
    112 ! the lower left corner of the grid is used.                             *
    113 ! Unlike in the eta system, no difference between heights for u,v and    *
    114 ! heights for w exists.                                                  *
     105! u,v,T and qv are given.                                                *
    115106!*************************************************************************
    116 
    117 
    118 !eso measure CPU time
    119 !  call mpif_mtime('verttransform',0)
    120107
    121108  if (init) then
    … …  
    123110! Search for a point with high surface pressure (i.e. not above significant topography)
    124111! Then, use this point to construct a reference z profile, to be used at all times
    125 !*****************************************************************************
     112!**************************************************************************************
    126113
    127114    do jy=0,nymin1
    … …  
    137124
    138125
    139     tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
     126    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    140127         ps(ixm,jym,1,n))
    141     pold(ixm,jym)=ps(ixm,jym,1,n)
     128    pold=ps(ixm,jym,1,n)
    142129    height(1)=0.
    143130
    … …  
    146133      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    147134
    148       if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
    149         height(kz)= &
    150              height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
    151              (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
     135      if (abs(tv-tvold).gt.0.2) then
     136        height(kz)=height(kz-1)+const*log(pold/pint)* &
     137             (tv-tvold)/log(tv/tvold)
    152138      else
    153         height(kz)=height(kz-1)+ &
    154              const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
     139        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    155140      endif
    156141
    157       tvold(ixm,jym)=tv(ixm,jym)
    158       pold(ixm,jym)=pint(ixm,jym)
     142      tvold=tv
     143      pold=pint
    159144    end do
    160145
    … …  
    182167!*************************
    183168
    184 
    185169  do jy=0,nymin1
    186170    do ix=0,nxmin1
    187       tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
    188            ps(ix,jy,1,n))
    189     enddo
    190   enddo
    191   pold=ps(:,:,1,n)
    192   uvzlev(:,:,1)=0.
    193   wzlev(:,:,1)=0.
    194   rhoh(:,:,1)=pold/(r_air*tvold)
     171      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
     172      pold=ps(ix,jy,1,n)
     173      uvwzlev(ix,jy,1)=0.
     174      wzlev(1)=0.
     175      rhoh(1)=pold/(r_air*tvold)
    195176
    196177
    … …  
    198179!******************************
    199180
    200   do kz=2,nuvz
    201     pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
    202     tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
    203     rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
    204 
    205     where (abs(tv-tvold).gt.0.2)
    206       uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
    207            (tv-tvold)/log(tv/tvold)
    208     elsewhere
    209       uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
    210     endwhere
    211 
    212     tvold=tv
    213     pold=pint
    214   end do
    215 
    216 
    217   do kz=2,nwz-1
    218     wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
    219   end do
    220   wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
    221        uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
     181      do kz=2,nuvz
     182        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
     183        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
     184        rhoh(kz)=pint/(r_air*tv)
     185
     186        if (abs(tv-tvold).gt.0.2) then
     187          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     188               (tv-tvold)/log(tv/tvold)
     189        else
     190          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
     191        endif
     192
     193        tvold=tv
     194        pold=pint
     195      end do
     196
     197
     198      do kz=2,nwz-1
     199        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
     200      end do
     201      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    222202
    223203! pinmconv=(h2-h1)/(p2-p1)
    224204
    225   pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
    226        ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
    227        (aknew(1)+bknew(1)*ps(:,:,1,n)))
    228   do kz=2,nz-1
    229     pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
    230          ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
    231          (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
    232   end do
    233   pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
    234        ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
    235        (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
     205      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
     206           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     207           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     208      do kz=2,nz-1
     209        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
     210             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     211             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     212      end do
     213      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
     214           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     215           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    236216
    237217! Levels, where u,v,t and q are given
    238218!************************************
    239219
    240 
    241   uu(:,:,1,n)=uuh(:,:,1)
    242   vv(:,:,1,n)=vvh(:,:,1)
    243   tt(:,:,1,n)=tth(:,:,1,n)
    244   qv(:,:,1,n)=qvh(:,:,1,n)
     220      uu(ix,jy,1,n)=uuh(ix,jy,1)
     221      vv(ix,jy,1,n)=vvh(ix,jy,1)
     222      tt(ix,jy,1,n)=tth(ix,jy,1,n)
     223      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
    245224!hg adding the cloud water
    246   clwc(:,:,1,n)=clwch(:,:,1,n)
    247   ciwc(:,:,1,n)=ciwch(:,:,1,n)   
     225      clwc(ix,jy,1,n)=clwch(ix,jy,1,n)
     226      ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)   
    248227!hg
    249   pv(:,:,1,n)=pvh(:,:,1)
    250   rho(:,:,1,n)=rhoh(:,:,1)
    251   uu(:,:,nz,n)=uuh(:,:,nuvz)
    252   vv(:,:,nz,n)=vvh(:,:,nuvz)
    253   tt(:,:,nz,n)=tth(:,:,nuvz,n)
    254   qv(:,:,nz,n)=qvh(:,:,nuvz,n)
     228      pv(ix,jy,1,n)=pvh(ix,jy,1)
     229      rho(ix,jy,1,n)=rhoh(1)
     230      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
     231      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
     232      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
     233      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
    255234
    256235!hg adding the cloud water
    257   clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
    258   ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
     236      clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
     237      ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)
    259238!hg
    260   pv(:,:,nz,n)=pvh(:,:,nuvz)
    261   rho(:,:,nz,n)=rhoh(:,:,nuvz)
    262 
    263 
    264   kmin=2
    265   idx=kmin
    266   do iz=2,nz-1
    267     do jy=0,nymin1
    268       do ix=0,nxmin1
    269         if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
    270           uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    271           vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    272           tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
    273           qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
     239      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
     240      rho(ix,jy,nz,n)=rhoh(nuvz)
     241      kmin=2
     242      do iz=2,nz-1
     243        do kz=kmin,nuvz
     244          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     245            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
     246            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
     247            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
     248            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
    274249!hg adding the cloud water
    275           clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
    276           ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     250            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
     251            ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
    277252!hg
    278           pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
    279           rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
    280         else
    281           innuvz: do kz=idx(ix,jy),nuvz
    282             if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    283                  (height(iz).le.uvzlev(ix,jy,kz))) then
    284               idx(ix,jy)=kz
    285               exit innuvz
    286             endif
    287           enddo innuvz
    288         endif
    289       enddo
    290     enddo
    291     do jy=0,nymin1
    292       do ix=0,nxmin1
    293         if(height(iz).le.uvzlev(ix,jy,nuvz)) then
    294           kz=idx(ix,jy)
    295           dz1=height(iz)-uvzlev(ix,jy,kz-1)
    296           dz2=uvzlev(ix,jy,kz)-height(iz)
    297           dz=dz1+dz2
    298           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    299           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
    300           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
    301                +tth(ix,jy,kz,n)*dz1)/dz
    302           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
    303                +qvh(ix,jy,kz,n)*dz1)/dz
     253            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
     254            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
     255            goto 30
     256          endif
     257          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     258               (height(iz).le.uvwzlev(ix,jy,kz))) then
     259            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     260            dz2=uvwzlev(ix,jy,kz)-height(iz)
     261            dz=dz1+dz2
     262            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
     263            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
     264            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2+tth(ix,jy,kz,n)*dz1)/dz
     265            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2+qvh(ix,jy,kz,n)*dz1)/dz
    304266!hg adding the cloud water
    305           clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
    306           ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
     267            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
     268            ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
    307269!hg
    308           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    309           rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
    310         endif
    311       enddo
    312     enddo
    313   enddo
     270
     271            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
     272            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
     273            kmin=kz
     274            goto 30
     275          endif
     276        end do
     27730      continue
     278      end do
    314279
    315280
    … …  
    317282!*************************
    318283
    319   ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
    320   ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
    321   kmin=2
    322   idx=kmin
    323   do iz=2,nz
    324     do jy=0,nymin1
    325       do ix=0,nxmin1
    326         inn:         do kz=idx(ix,jy),nwz
    327           if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
    328                height(iz).le.wzlev(ix,jy,kz)) then
    329             idx(ix,jy)=kz
    330             exit inn
     284      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
     285      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
     286      kmin=2
     287      do iz=2,nz
     288        do kz=kmin,nwz
     289          if ((height(iz).gt.wzlev(kz-1)).and. &
     290               (height(iz).le.wzlev(kz))) then
     291            dz1=height(iz)-wzlev(kz-1)
     292            dz2=wzlev(kz)-height(iz)
     293            dz=dz1+dz2
     294            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
     295                 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
     296            kmin=kz
     297            goto 40
    331298          endif
    332         enddo inn
    333       enddo
    334     enddo
    335     do jy=0,nymin1
    336       do ix=0,nxmin1
    337         kz=idx(ix,jy)
    338         dz1=height(iz)-wzlev(ix,jy,kz-1)
    339         dz2=wzlev(ix,jy,kz)-height(iz)
    340         dz=dz1+dz2
    341         ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
    342              +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
    343       enddo
    344     enddo
    345   enddo
     299        end do
     30040      continue
     301      end do
    346302
    347303! Compute density gradients at intermediate levels
    348304!*************************************************
    349305
    350   drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
    351        (height(2)-height(1))
    352   do kz=2,nz-1
    353     drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
    354          (height(kz+1)-height(kz-1))
     306      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
     307           (height(2)-height(1))
     308      do kz=2,nz-1
     309        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
     310             (height(kz+1)-height(kz-1))
     311      end do
     312      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
     313
     314    end do
    355315  end do
    356   drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
    357 
    358 !    end do
    359 !  end do
    360316
    361317
    … …  
    366322
    367323  do jy=1,ny-2
    368     cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
    369   enddo
    370 
    371   kmin=2
    372   idx=kmin
    373   do iz=2,nz-1
    374     do jy=1,ny-2
    375       do ix=1,nx-2
    376 
    377         inneta: do kz=idx(ix,jy),nz
    378           if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
    379                (height(iz).le.uvzlev(ix,jy,kz))) then
    380             idx(ix,jy)=kz
    381             exit inneta
     324    cosf=cos((real(jy)*dy+ylat0)*pi180)
     325    do ix=1,nx-2
     326
     327      kmin=2
     328      do iz=2,nz-1
     329
     330        ui=uu(ix,jy,iz,n)*dxconst/cosf
     331        vi=vv(ix,jy,iz,n)*dyconst
     332
     333        do kz=kmin,nz
     334          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     335               (height(iz).le.uvwzlev(ix,jy,kz))) then
     336            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     337            dz2=uvwzlev(ix,jy,kz)-height(iz)
     338            dz=dz1+dz2
     339            kl=kz-1
     340            klp=kz
     341            kmin=kz
     342            goto 47
    382343          endif
    383         enddo inneta
    384       enddo
    385     enddo
    386 
    387     do jy=1,ny-2
    388       do ix=1,nx-2
    389         kz=idx(ix,jy)
    390         dz1=height(iz)-uvzlev(ix,jy,kz-1)
    391         dz2=uvzlev(ix,jy,kz)-height(iz)
    392         dz=dz1+dz2
    393         ix1=ix-1
     344        end do
     345
     34647      ix1=ix-1
    394347        jy1=jy-1
    395348        ixp=ix+1
    396349        jyp=jy+1
    397350
    398         dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
    399         dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
     351        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
     352        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
    400353        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
    401354
    402         dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
    403         dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
     355        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
     356        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
    404357        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    405358
    406         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)
     359        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
    407360
    408361      end do
    … …  
    410363    end do
    411364  end do
     365
    412366
    413367! If north pole is in the domain, calculate wind velocities in polar
    … …  
    416370
    417371  if (nglobal) then
    418     do iz=1,nz
    419       do jy=int(switchnorthg)-2,nymin1
    420         ylat=ylat0+real(jy)*dy
    421         do ix=0,nxmin1
    422           xlon=xlon0+real(ix)*dx
     372    do jy=int(switchnorthg)-2,nymin1
     373      ylat=ylat0+real(jy)*dy
     374      do ix=0,nxmin1
     375        xlon=xlon0+real(ix)*dx
     376        do iz=1,nz
    423377          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    424                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    425                vvpol(ix,jy,iz,n))
     378               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    426379        end do
    427380      end do
    … …  
    437390      xlon=xlon0+real(nx/2-1)*dx
    438391      xlonr=xlon*pi/180.
    439       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
    440            vv(nx/2-1,nymin1,iz,n)**2)
     392      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
    441393      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
    442         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
    443              vv(nx/2-1,nymin1,iz,n))-xlonr
     394        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
    444395      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    445396        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    … …  
    457408      uuaux=-ffpol*sin(xlonr+ddpol)
    458409      vvaux=-ffpol*cos(xlonr+ddpol)
    459       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    460            vvpolaux)
    461 
     410      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    462411      jy=nymin1
    463412      do ix=0,nxmin1
    … …  
    492441
    493442  if (sglobal) then
    494     do iz=1,nz
    495       do jy=0,int(switchsouthg)+3
    496         ylat=ylat0+real(jy)*dy
    497         do ix=0,nxmin1
    498           xlon=xlon0+real(ix)*dx
     443    do jy=0,int(switchsouthg)+3
     444      ylat=ylat0+real(jy)*dy
     445      do ix=0,nxmin1
     446        xlon=xlon0+real(ix)*dx
     447        do iz=1,nz
    499448          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    500                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    501                vvpol(ix,jy,iz,n))
     449               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    502450        end do
    503451      end do
    … …  
    512460      xlon=xlon0+real(nx/2-1)*dx
    513461      xlonr=xlon*pi/180.
    514       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
    515            vv(nx/2-1,0,iz,n)**2)
     462      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
    516463      if (vv(nx/2-1,0,iz,n).lt.0.) then
    517         ddpol=atan(uu(nx/2-1,0,iz,n)/ &
    518              vv(nx/2-1,0,iz,n))+xlonr
     464        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    519465      else if (vv(nx/2-1,0,iz,n).gt.0.) then
    520         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
    521              vv(nx/2-1,0,iz,n))+xlonr
     466        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    522467      else
    523468        ddpol=pi/2-xlonr
    … …  
    532477      uuaux=+ffpol*sin(xlonr-ddpol)
    533478      vvaux=-ffpol*cos(xlonr-ddpol)
    534       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    535            vvpolaux)
     479      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    536480
    537481      jy=0
    … …  
    561505
    562506
     507!write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
     508!   create a cloud and rainout/washout field, clouds occur where rh>80%
     509!   total cloudheight is stored at level 0
     510
    563511  if (readclouds) write(*,*) 'using cloud water from ECMWF'
    564512!  if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
    … …  
    570518    do ix=0,nxmin1
    571519! OLD METHOD
    572       rain_cloud_above(ix,jy)=0
     520      rain_cloud_above=0
    573521      lsp=lsprec(ix,jy,1,n)
    574522      convp=convprec(ix,jy,1,n)
    … …  
    581529        if (rh.gt.0.8) then ! in cloud
    582530          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    583             rain_cloud_above(ix,jy)=1
     531            rain_cloud_above=1
    584532            cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    585533                 height(kz)-height(kz-1)
    … …  
    593541          endif
    594542        else ! no cloud
    595           if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
     543          if (rain_cloud_above.eq.1) then ! scavenging
    596544            if (lsp.ge.convp) then
    597545              clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    … …  
    605553
    606554! PS 2012
    607 !      lsp=lsprec(ix,jy,1,n)
    608 !      convp=convprec(ix,jy,1,n)
     555!          lsp=lsprec(ix,jy,1,n)
     556!          convp=convprec(ix,jy,1,n)
    609557!          prec=lsp+convp
    610558!          if (lsp.gt.convp) then !  prectype='lsp'
    … …  
    612560!          else ! prectype='cp '
    613561!            lconvectprec = .true.
    614 !           endif
     562!          endif
    615563!HG METHOD
    616564!readclouds =.true.
    … …  
    624572!        do kz=1, nz
    625573!         !clw & ciw are given in kg/kg 
    626           ! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
    627           ! if (cloud_water .gt. 0) then
    628           !   cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
    629           !   cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
    630           !   cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid
    631           ! endif
    632           ! cloud_ver=max(0,cloud_max-cloud_min)
    633 
    634           ! icloudbot(ix,jy,n)=cloud_min
    635           ! icloudthck(ix,jy,n)=cloud_ver
    636           ! rcw(ix,jy)=cloud_col_wat
    637           ! rpc(ix,jy)=prec
     574!         cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
     575!         if (cloud_water .gt. 0) then
     576!          cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
     577!          cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
     578!          cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid
     579!         endif
     580!         cloud_ver=max(0,cloud_max-cloud_min)
     581
     582! if (clwc(ix,jy,kz,n).gt.0 .or.  ciwc(ix,jy,kz,n).gt.0) &
     583!     !write(*,*) 'WATER',clwc(ix,jy,kz,n), 'ICE',ciwc(ix,jy,kz,n),'RH',rh,'KZ',kz &
     584!     write(*,*) 'WATER',cloud_water,'RH',rh,'PREC',prec,'HEIGHT',height(kz) &
     585!       enddo
     586!       if ( cloud_min .ne. 99999 .and. cloud_max .ne. -1) write(*,*) 'CB', cloud_min, '   CT',cloud_max
     587! if (prec .gt. 0) write(*,*) 'PREC',prec,'Cloud Bot',cloud_min,'Cloud Top',cloud_max, 'Cloud Vert. ext',cloud_ver &
     588!                             ,'COLUMN cloud water',cloud_col_wat,'Conevctive' ,lconvectprec
     589!       icloudbot(ix,jy,n)=cloud_min
     590!       icloudthck(ix,jy,n)=cloud_ver
     591!       rcw(ix,jy)=cloud_col_wat
     592!       rpc(ix,jy)=prec
    638593!write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code
    639594!END HG METHOD
    … …  
    658613!              icloudtop=nint(height(kz)) ! use int to save memory
    659614!            endif
    660     ! end do
     615!          enddo
    661616!PS try to get a cloud thicker than 50 m
    662617!PS if there is at least .01 mm/h  - changed to 0.002 and put into
    … …  
    667622!            rhmin = rhmin - 0.05
    668623!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
    669 !   end if
     624!          endif
    670625
    671626!PS is based on looking at a limited set of comparison data
    … …  
    701656  end do
    702657
    703 !do 102 kz=1,nuvz
    704 !write(an,'(i02)') kz+10
    705 !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
    706 !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
    707 !do 101 jy=0,nymin1
    708 !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
    709 !101   continue
    710 ! close(4)
    711 !102   continue
    712 
    713 ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
    714 ! do 103 jy=0,nymin1
    715 !     write (4,*)
    716 !+       (height(kz),kz=1,nuvz)
    717 !103   continue
    718 ! close(4)
    719 
    720 !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
    721 ! do 104 jy=0,nymin1
    722 !     write (4,*)
    723 !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
    724 !104   continue
    725 ! close(4)
    726 
    727 !eso measure CPU time
    728 !  call mpif_mtime('verttransform',1)
    729 
    730 !eso print out the same measure as in Leo's routine
    731     ! write(*,*) 'verttransform: ', &
    732     !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
    733     !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
    734     !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
    735     !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
    736     !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
    737     !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
    738   end subroutine verttransform
    739 
     658end subroutine verttransform
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG