Changeset 8a65cb0 in flexpart.git for src/verttransform.f90


Ignore:
Timestamp:
Mar 2, 2015, 3:11:55 PM (9 years ago)
Author:
Espen Sollum ATMOS <espen@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
1d207bb
Parents:
60403cd
Message:

Added code, makefile for dev branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/verttransform.f90

    r4fbe7a5 r8a65cb0  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    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   !                                                                            *
    33   !     Author: A. Stohl, G. Wotawa                                            *
    34   !                                                                            *
    35   !     12 August 1996                                                         *
    36   !     Update: 16 January 1998                                                *
    37   !                                                                            *
    38   !     Major update: 17 February 1999                                         *
    39   !     by G. Wotawa                                                           *
    40   !                                                                            *
    41   !     - Vertical levels for u, v and w are put together                      *
    42   !     - Slope correction for vertical velocity: Modification of calculation  *
    43   !       procedure                                                            *
    44   !                                                                            *
    45   !*****************************************************************************
    46   !  Changes, Bernd C. Krueger, Feb. 2001:
    47   !   Variables tth and qvh (on eta coordinates) from common block
    48   !*****************************************************************************
    49   ! Sabine Eckhardt, March 2007
    50   ! added the variable cloud for use with scavenging - descr. in com_mod
    51   !*****************************************************************************
    52   !                                                                            *
    53   ! Variables:                                                                 *
    54   ! nx,ny,nz                        field dimensions in x,y and z direction    *
    55   ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
    56   ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
    57   ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
    58   ! ww(0:nxmax,0:nymax,nzmax,2)     wind components in z-direction [deltaeta/s]*
    59   ! tt(0:nxmax,0:nymax,nzmax,2)     temperature [K]                            *
    60   ! pv(0:nxmax,0:nymax,nzmax,2)     potential voriticity (pvu)                 *
    61   ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
    62   !                                                                            *
    63   !*****************************************************************************
     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!                                                                            *
     33!     Author: A. Stohl, G. Wotawa                                            *
     34!                                                                            *
     35!     12 August 1996                                                         *
     36!     Update: 16 January 1998                                                *
     37!                                                                            *
     38!     Major update: 17 February 1999                                         *
     39!     by G. Wotawa                                                           *
     40!                                                                            *
     41!     - Vertical levels for u, v and w are put together                      *
     42!     - Slope correction for vertical velocity: Modification of calculation  *
     43!       procedure                                                            *
     44!                                                                            *
     45!*****************************************************************************
     46!  Changes, Bernd C. Krueger, Feb. 2001:
     47!   Variables tth and qvh (on eta coordinates) from common block
     48!*****************************************************************************
     49! Sabine Eckhardt, March 2007
     50! added the variable cloud for use with scavenging - descr. in com_mod
     51!*****************************************************************************
     52!                                                                            *
     53! Variables:                                                                 *
     54! nx,ny,nz                        field dimensions in x,y and z direction    *
     55! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition    *
     56! uu(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in x-direction [m/s]*
     57! vv(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in y-direction [m/s]*
     58! ww(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in z-direction      *
     59!                                          [deltaeta/s]                      *
     60! tt(0:nxmax,0:nymax,nzmax,numwfmem)     temperature [K]                     *
     61! pv(0:nxmax,0:nymax,nzmax,numwfmem)     potential voriticity (pvu)          *
     62! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
     63!                                                                            *
     64!*****************************************************************************
    6465
    6566  use par_mod
    … …  
    7172  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    7273  integer :: rain_cloud_above,kz_inv
     74! integer :: icloudtop !PS
    7375  real :: f_qvsat,pressure
    74   real :: rh,lsp,convp
     76! real :: rh,lsp,convp
     77  real :: rh,lsp,convp,prec,rhmin,precmin
    7578  real :: rhoh(nuvzmax),pinmconv(nzmax)
    7679  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    … …  
    8386  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8487  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
     88  logical lconvectprec
    8589  real,parameter :: const=r_air/ga
     90  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    8691
    8792  logical :: init = .true.
    8893
    89 
    90   !*************************************************************************
    91   ! If verttransform is called the first time, initialize heights of the   *
    92   ! z levels in meter. The heights are the heights of model levels, where  *
    93   ! u,v,T and qv are given.                                                *
    94   !*************************************************************************
     94!hg
     95  integer :: cloud_ver,cloud_min, cloud_max
     96  real :: cloud_col_wat, cloud_water
     97!hg temporary variables for testing
     98  real :: rcw(0:nxmax-1,0:nymax-1)
     99  real :: rpc(0:nxmax-1,0:nymax-1)
     100!hg
     101
     102!*************************************************************************
     103! If verttransform is called the first time, initialize heights of the   *
     104! z levels in meter. The heights are the heights of model levels, where  *
     105! u,v,T and qv are given.                                                *
     106!*************************************************************************
    95107
    96108  if (init) then
    97109
    98   ! Search for a point with high surface pressure (i.e. not above significant topography)
    99   ! Then, use this point to construct a reference z profile, to be used at all times
    100   !**************************************************************************************
     110! Search for a point with high surface pressure (i.e. not above significant topography)
     111! Then, use this point to construct a reference z profile, to be used at all times
     112!**************************************************************************************
    101113
    102114    do jy=0,nymin1
    … …  
    113125
    114126    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    115     ps(ixm,jym,1,n))
     127         ps(ixm,jym,1,n))
    116128    pold=ps(ixm,jym,1,n)
    117129    height(1)=0.
    … …  
    123135      if (abs(tv-tvold).gt.0.2) then
    124136        height(kz)=height(kz-1)+const*log(pold/pint)* &
    125         (tv-tvold)/log(tv/tvold)
     137             (tv-tvold)/log(tv/tvold)
    126138      else
    127139        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    … …  
    133145
    134146
    135   ! Determine highest levels that can be within PBL
    136   !************************************************
     147! Determine highest levels that can be within PBL
     148!************************************************
    137149
    138150    do kz=1,nz
    … …  
    1441569   continue
    145157
    146   ! Do not repeat initialization of the Cartesian z grid
    147   !*****************************************************
     158! Do not repeat initialization of the Cartesian z grid
     159!*****************************************************
    148160
    149161    init=.false.
    … …  
    152164
    153165
    154   ! Loop over the whole grid
    155   !*************************
     166! Loop over the whole grid
     167!*************************
    156168
    157169  do jy=0,nymin1
    … …  
    164176
    165177
    166   ! Compute heights of eta levels
    167   !******************************
     178! Compute heights of eta levels
     179!******************************
    168180
    169181      do kz=2,nuvz
    … …  
    174186        if (abs(tv-tvold).gt.0.2) then
    175187          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    176           (tv-tvold)/log(tv/tvold)
     188               (tv-tvold)/log(tv/tvold)
    177189        else
    178190          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    … …  
    189201      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    190202
    191   ! pinmconv=(h2-h1)/(p2-p1)
     203! pinmconv=(h2-h1)/(p2-p1)
    192204
    193205      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    194       ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
    195       (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     206           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     207           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
    196208      do kz=2,nz-1
    197209        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    198         ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
    199         (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     210             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     211             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
    200212      end do
    201213      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    202       ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
    203       (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    204 
    205   ! Levels, where u,v,t and q are given
    206   !************************************
     214           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     215           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
     216
     217! Levels, where u,v,t and q are given
     218!************************************
    207219
    208220      uu(ix,jy,1,n)=uuh(ix,jy,1)
    … …  
    210222      tt(ix,jy,1,n)=tth(ix,jy,1,n)
    211223      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
     224!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)   
     227!hg
    212228      pv(ix,jy,1,n)=pvh(ix,jy,1)
    213229      rho(ix,jy,1,n)=rhoh(1)
    … …  
    216232      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
    217233      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
     234
     235!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)
     238!hg
    218239      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
    219240      rho(ix,jy,nz,n)=rhoh(nuvz)
    … …  
    226247            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
    227248            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
     249!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)
     252!hg
    228253            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
    229254            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
    … …  
    231256          endif
    232257          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    233           (height(iz).le.uvwzlev(ix,jy,kz))) then
    234            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    235            dz2=uvwzlev(ix,jy,kz)-height(iz)
    236            dz=dz1+dz2
    237            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    238            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
    239            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
    240                 +tth(ix,jy,kz,n)*dz1)/dz
    241            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
    242                 +qvh(ix,jy,kz,n)*dz1)/dz
    243            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    244            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    245            kmin=kz
    246            goto 30
     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
     266!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
     269!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
    247275          endif
    248276        end do
    … …  
    251279
    252280
    253   ! Levels, where w is given
    254   !*************************
     281! Levels, where w is given
     282!*************************
    255283
    256284      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
    … …  
    261289          if ((height(iz).gt.wzlev(kz-1)).and. &
    262290               (height(iz).le.wzlev(kz))) then
    263            dz1=height(iz)-wzlev(kz-1)
    264            dz2=wzlev(kz)-height(iz)
    265            dz=dz1+dz2
    266            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
    267                 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
    268            kmin=kz
    269            goto 40
     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
    270298          endif
    271299        end do
    … …  
    273301      end do
    274302
    275   ! Compute density gradients at intermediate levels
    276   !*************************************************
     303! Compute density gradients at intermediate levels
     304!*************************************************
    277305
    278306      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
    279       (height(2)-height(1))
     307           (height(2)-height(1))
    280308      do kz=2,nz-1
    281309        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    282         (height(kz+1)-height(kz-1))
     310             (height(kz+1)-height(kz-1))
    283311      end do
    284312      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    288316
    289317
    290   !****************************************************************
    291   ! Compute slope of eta levels in windward direction and resulting
    292   ! vertical wind correction
    293   !****************************************************************
     318!****************************************************************
     319! Compute slope of eta levels in windward direction and resulting
     320! vertical wind correction
     321!****************************************************************
    294322
    295323  do jy=1,ny-2
    … …  
    305333        do kz=kmin,nz
    306334          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    307           (height(iz).le.uvwzlev(ix,jy,kz))) then
     335               (height(iz).le.uvwzlev(ix,jy,kz))) then
    308336            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    309337            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    337365
    338366
    339   ! If north pole is in the domain, calculate wind velocities in polar
    340   ! stereographic coordinates
    341   !*******************************************************************
     367! If north pole is in the domain, calculate wind velocities in polar
     368! stereographic coordinates
     369!*******************************************************************
    342370
    343371  if (nglobal) then
    … …  
    348376        do iz=1,nz
    349377          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    350           vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     378               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    351379        end do
    352380      end do
    … …  
    356384    do iz=1,nz
    357385
    358   ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
    359   !
    360   !   AMSnauffer Nov 18 2004 Added check for case vv=0
    361   !
     386! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
     387!
     388!   AMSnauffer Nov 18 2004 Added check for case vv=0
     389!
    362390      xlon=xlon0+real(nx/2-1)*dx
    363391      xlonr=xlon*pi/180.
    … …  
    367395      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    368396        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    369         vv(nx/2-1,nymin1,iz,n))-xlonr
     397             vv(nx/2-1,nymin1,iz,n))-xlonr
    370398      else
    371399        ddpol=pi/2-xlonr
    … …  
    374402      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
    375403
    376   ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
     404! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
    377405      xlon=180.0
    378406      xlonr=xlon*pi/180.
    … …  
    389417
    390418
    391   ! Fix: Set W at pole to the zonally averaged W of the next equator-
    392   ! ward parallel of latitude
    393 
    394   do iz=1,nz
     419! Fix: Set W at pole to the zonally averaged W of the next equator-
     420! ward parallel of latitude
     421
     422    do iz=1,nz
    395423      wdummy=0.
    396424      jy=ny-2
    … …  
    403431        ww(ix,jy,iz,n)=wdummy
    404432      end do
    405   end do
     433    end do
    406434
    407435  endif
    408436
    409437
    410   ! If south pole is in the domain, calculate wind velocities in polar
    411   ! stereographic coordinates
    412   !*******************************************************************
     438! If south pole is in the domain, calculate wind velocities in polar
     439! stereographic coordinates
     440!*******************************************************************
    413441
    414442  if (sglobal) then
    … …  
    419447        do iz=1,nz
    420448          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    421           vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     449               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    422450        end do
    423451      end do
    … …  
    426454    do iz=1,nz
    427455
    428   ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
    429   !
    430   !   AMSnauffer Nov 18 2004 Added check for case vv=0
    431   !
     456! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
     457!
     458!   AMSnauffer Nov 18 2004 Added check for case vv=0
     459!
    432460      xlon=xlon0+real(nx/2-1)*dx
    433461      xlonr=xlon*pi/180.
    … …  
    443471      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
    444472
    445   ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
     473! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
    446474      xlon=180.0
    447475      xlonr=xlon*pi/180.
    … …  
    459487
    460488
    461   ! Fix: Set W at pole to the zonally averaged W of the next equator-
    462   ! ward parallel of latitude
     489! Fix: Set W at pole to the zonally averaged W of the next equator-
     490! ward parallel of latitude
    463491
    464492    do iz=1,nz
    … …  
    477505
    478506
    479   !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
    480   !   create a cloud and rainout/washout field, clouds occur where rh>80%
    481   !   total cloudheight is stored at level 0
     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
     511  if (readclouds) write(*,*) 'using cloud water from ECMWF'
     512  if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
     513
     514  rcw(:,:)=0
     515  rpc(:,:)=0
     516
    482517  do jy=0,nymin1
    483518    do ix=0,nxmin1
     519! OLD METHOD
    484520      rain_cloud_above=0
    485521      lsp=lsprec(ix,jy,1,n)
    … …  
    487523      cloudsh(ix,jy,n)=0
    488524      do kz_inv=1,nz-1
    489          kz=nz-kz_inv+1
    490          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    491          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    492          clouds(ix,jy,kz,n)=0
    493          if (rh.gt.0.8) then ! in cloud
    494             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    495                rain_cloud_above=1
    496                cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    497                height(kz)-height(kz-1)
    498                if (lsp.ge.convp) then
    499                   clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    500                else
    501                   clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    502                endif
    503             else ! no precipitation
    504                   clouds(ix,jy,kz,n)=1 ! cloud
     525        kz=nz-kz_inv+1
     526        pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     527        rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     528        clouds(ix,jy,kz,n)=0
     529        if (rh.gt.0.8) then ! in cloud
     530          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     531            rain_cloud_above=1
     532            cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
     533                 height(kz)-height(kz-1)
     534            if (lsp.ge.convp) then
     535              clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
     536            else
     537              clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    505538            endif
    506          else ! no cloud
    507             if (rain_cloud_above.eq.1) then ! scavenging
    508                if (lsp.ge.convp) then
    509                   clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    510                else
    511                   clouds(ix,jy,kz,n)=4 ! convp dominated washout
    512                endif
     539          else ! no precipitation
     540            clouds(ix,jy,kz,n)=1 ! cloud
     541          endif
     542        else ! no cloud
     543          if (rain_cloud_above.eq.1) then ! scavenging
     544            if (lsp.ge.convp) then
     545              clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     546            else
     547              clouds(ix,jy,kz,n)=4 ! convp dominated washout
    513548            endif
    514          endif
    515       end do
     549          endif
     550        endif
     551      end do
     552!END OLD METHOD
     553
     554! PS 2012
     555!          lsp=lsprec(ix,jy,1,n)
     556!          convp=convprec(ix,jy,1,n)
     557!          prec=lsp+convp
     558!          if (lsp.gt.convp) then !  prectype='lsp'
     559!            lconvectprec = .false.
     560!          else ! prectype='cp '
     561!            lconvectprec = .true.
     562!          endif
     563!HG METHOD
     564!readclouds =.true.
     565!      if (readclouds) then
     566!hg added APR 2014  Cloud Water=clwc(ix,jy,kz,n)  Cloud Ice=ciwc(ix,jy,kz,n)
     567!hg Use the cloud water variables to determine existence of clouds. This makes the PS code obsolete
     568!        cloud_min=99999
     569!        cloud_max=-1
     570!        cloud_col_wat=0
     571
     572!        do kz=1, nz
     573!         !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
     593!write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code
     594!END HG METHOD
     595
     596
     597
     598!      else ! windfields does not contain cloud data
     599!          rhmin = 0.90 ! standard condition for presence of clouds
     600!PS       note that original by Sabine Eckhart was 80%
     601!PS       however, for T<-20 C we consider saturation over ice
     602!PS       so I think 90% should be enough         
     603!          icloudbot(ix,jy,n)=icmv
     604!          icloudtop=icmv ! this is just a local variable
     605!98        do kz=1,nz
     606!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     607!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     608!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
     609!            if (rh .gt. rhmin) then
     610!              if (icloudbot(ix,jy,n) .eq. icmv) then
     611!                icloudbot(ix,jy,n)=nint(height(kz))
     612!              endif
     613!              icloudtop=nint(height(kz)) ! use int to save memory
     614!            endif
     615!          enddo
     616!PS try to get a cloud thicker than 50 m
     617!PS if there is at least .01 mm/h  - changed to 0.002 and put into
     618!PS parameter precpmin       
     619!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
     620!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
     621!              prec .gt. precmin) then
     622!            rhmin = rhmin - 0.05
     623!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
     624!          endif
     625
     626!PS is based on looking at a limited set of comparison data
     627!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
     628!             prec .gt. precmin) then
     629!
     630!            if (convp .lt. 0.1) then
     631!              icloudbot(ix,jy,n) = 500
     632!              icloudtop =         8000
     633!            else
     634!              icloudbot(ix,jy,n) = 0
     635!              icloudtop =      10000
     636!            endif
     637!          endif
     638!          if (icloudtop .ne. icmv) then
     639!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
     640!          else
     641!            icloudthck(ix,jy,n) = icmv
     642!          endif
     643!PS  get rid of too thin clouds     
     644!          if (icloudthck(ix,jy,n) .lt. 50) then
     645!            icloudbot(ix,jy,n)=icmv
     646!            icloudthck(ix,jy,n)=icmv
     647!          endif
     648!hg__________________________________
     649!           rcw(ix,jy)=2E-7*prec**0.36
     650!           rpc(ix,jy)=prec
     651!hg end______________________________
     652
     653!      endif !hg read clouds
     654
    516655    end do
    517656  end do
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG