Changeset 4d45639 in flexpart.git


Ignore:
Timestamp:
Nov 20, 2015, 4:00:21 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:
b255cd0
Parents:
adf46ae
Message:

Merged in optimized verttransform.f90 from author Leopold Haimberger

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/verttransform.f90

    r0f20c31 r4d45639  
    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
    6974
    7075  implicit none
    7176
    7277  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    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
     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
    8085  real :: xlon,ylat,xlonr,dzdx,dzdy
    81   real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
     86  real :: dzdx1,dzdx2,dzdy1,dzdy2
    8287  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8388  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     
    8590  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
    8691  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    87   real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
    88   logical lconvectprec
     92  real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)
    8993  real,parameter :: const=r_air/ga
     94!  integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2
    9095  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9196
     
    103108! If verttransform is called the first time, initialize heights of the   *
    104109! z levels in meter. The heights are the heights of model levels, where  *
    105 ! u,v,T and qv are given.                                                *
     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.                                                  *
    106115!*************************************************************************
     116
     117
     118!eso measure CPU time
     119!  call mpif_mtime('verttransform',0)
    107120
    108121  if (init) then
     
    110123! Search for a point with high surface pressure (i.e. not above significant topography)
    111124! Then, use this point to construct a reference z profile, to be used at all times
    112 !**************************************************************************************
     125!*****************************************************************************
    113126
    114127    do jy=0,nymin1
     
    124137
    125138
    126     tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
     139    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    127140         ps(ixm,jym,1,n))
    128     pold=ps(ixm,jym,1,n)
     141    pold(ixm,jym)=ps(ixm,jym,1,n)
    129142    height(1)=0.
    130143
     
    133146      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    134147
    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)
     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))
    138152      else
    139         height(kz)=height(kz-1)+const*log(pold/pint)*tv
     153        height(kz)=height(kz-1)+ &
     154             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
    140155      endif
    141156
    142       tvold=tv
    143       pold=pint
     157      tvold(ixm,jym)=tv(ixm,jym)
     158      pold(ixm,jym)=pint(ixm,jym)
    144159    end do
    145160
     
    167182!*************************
    168183
     184
    169185  do jy=0,nymin1
    170186    do ix=0,nxmin1
    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)
     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)
    176195
    177196
     
    179198!******************************
    180199
    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)
     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)
    202222
    203223! pinmconv=(h2-h1)/(p2-p1)
    204224
    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)))
     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)))
    216236
    217237! Levels, where u,v,t and q are given
    218238!************************************
    219239
    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)
     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)
    224245!hg adding the cloud water
    225       clwc(ix,jy,1,n)=clwch(ix,jy,1,n)
    226       ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)   
     246  clwc(:,:,1,n)=clwch(:,:,1,n)
     247  ciwc(:,:,1,n)=ciwch(:,:,1,n)   
    227248!hg
    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)
     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)
    234255
    235256!hg adding the cloud water
    236       clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
    237       ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)
     257  clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
     258  ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
    238259!hg
    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)
     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)
    249274!hg adding the cloud water
    250             clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
    251             ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     275          clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
     276          ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
    252277!hg
    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
     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
    266304!hg adding the cloud water
    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
     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
    269307!hg
    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
    277 30      continue
    278       end do
     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
    279314
    280315
     
    282317!*************************
    283318
    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
     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
    298331          endif
    299         end do
    300 40      continue
    301       end do
     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
    302346
    303347! Compute density gradients at intermediate levels
    304348!*************************************************
    305349
    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
     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))
    315355  end do
     356  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
     357
     358!    end do
     359!  end do
    316360
    317361
     
    322366
    323367  do jy=1,ny-2
    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
     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
    343382          endif
    344         end do
    345 
    346 47      ix1=ix-1
     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
    347394        jy1=jy-1
    348395        ixp=ix+1
    349396        jyp=jy+1
    350397
    351         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
    352         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
     398        dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
     399        dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
    353400        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
    354401
    355         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
    356         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
     402        dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
     403        dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
    357404        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    358405
    359         ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
     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)
    360407
    361408      end do
     
    363410    end do
    364411  end do
    365 
    366412
    367413! If north pole is in the domain, calculate wind velocities in polar
     
    370416
    371417  if (nglobal) then
    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
     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
    377423          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    378                vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     424               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
     425               vvpol(ix,jy,iz,n))
    379426        end do
    380427      end do
     
    390437      xlon=xlon0+real(nx/2-1)*dx
    391438      xlonr=xlon*pi/180.
    392       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
     439      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
     440           vv(nx/2-1,nymin1,iz,n)**2)
    393441      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
    394         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
     442        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
     443             vv(nx/2-1,nymin1,iz,n))-xlonr
    395444      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    396445        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
     
    408457      uuaux=-ffpol*sin(xlonr+ddpol)
    409458      vvaux=-ffpol*cos(xlonr+ddpol)
    410       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     459      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     460           vvpolaux)
     461
    411462      jy=nymin1
    412463      do ix=0,nxmin1
     
    441492
    442493  if (sglobal) then
    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
     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
    448499          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    449                vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     500               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
     501               vvpol(ix,jy,iz,n))
    450502        end do
    451503      end do
     
    460512      xlon=xlon0+real(nx/2-1)*dx
    461513      xlonr=xlon*pi/180.
    462       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
     514      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
     515           vv(nx/2-1,0,iz,n)**2)
    463516      if (vv(nx/2-1,0,iz,n).lt.0.) then
    464         ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
     517        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
     518             vv(nx/2-1,0,iz,n))+xlonr
    465519      else if (vv(nx/2-1,0,iz,n).gt.0.) then
    466         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
     520        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
     521             vv(nx/2-1,0,iz,n))+xlonr
    467522      else
    468523        ddpol=pi/2-xlonr
     
    477532      uuaux=+ffpol*sin(xlonr-ddpol)
    478533      vvaux=-ffpol*cos(xlonr-ddpol)
    479       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     534      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     535           vvpolaux)
    480536
    481537      jy=0
     
    505561
    506562
    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 
    511563  if (readclouds) write(*,*) 'using cloud water from ECMWF'
    512564!  if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
     
    518570    do ix=0,nxmin1
    519571! OLD METHOD
    520       rain_cloud_above=0
     572      rain_cloud_above(ix,jy)=0
    521573      lsp=lsprec(ix,jy,1,n)
    522574      convp=convprec(ix,jy,1,n)
     
    529581        if (rh.gt.0.8) then ! in cloud
    530582          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    531             rain_cloud_above=1
     583            rain_cloud_above(ix,jy)=1
    532584            cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    533585                 height(kz)-height(kz-1)
     
    541593          endif
    542594        else ! no cloud
    543           if (rain_cloud_above.eq.1) then ! scavenging
     595          if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
    544596            if (lsp.ge.convp) then
    545597              clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     
    553605
    554606! PS 2012
    555 !          lsp=lsprec(ix,jy,1,n)
    556 !          convp=convprec(ix,jy,1,n)
     607!      lsp=lsprec(ix,jy,1,n)
     608!      convp=convprec(ix,jy,1,n)
    557609!          prec=lsp+convp
    558610!          if (lsp.gt.convp) then !  prectype='lsp'
     
    560612!          else ! prectype='cp '
    561613!            lconvectprec = .true.
    562 !          endif
     614!           endif
    563615!HG METHOD
    564616!readclouds =.true.
     
    572624!        do kz=1, nz
    573625!         !clw & ciw are given in kg/kg 
    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
     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
    593638!write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code
    594639!END HG METHOD
     
    613658!              icloudtop=nint(height(kz)) ! use int to save memory
    614659!            endif
    615 !          enddo
     660    ! end do
    616661!PS try to get a cloud thicker than 50 m
    617662!PS if there is at least .01 mm/h  - changed to 0.002 and put into
     
    622667!            rhmin = rhmin - 0.05
    623668!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
    624 !          endif
     669!   end if
    625670
    626671!PS is based on looking at a limited set of comparison data
     
    656701  end do
    657702
    658 end subroutine verttransform
     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
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG