pkg/atm_phys/lscale_cond_mod.F90

b2ea1d2979 Jean*0001 module lscale_cond_mod
                0002
                0003 !-----------------------------------------------------------------------
                0004 !use            fms_mod, only:  file_exist, error_mesg, open_file,  &
                0005 !                               check_nml_error, mpp_pe, FATAL,  &
                0006 !                               close_file
                0007 use simple_sat_vapor_pres_mod, only:  escomp, descomp
                0008
                0009 use      constants_mod, only:  HLv,HLs,Cp_air,Grav,rdgas,rvgas
                0010 use     gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0011
                0012 implicit none
                0013 private
                0014 !-----------------------------------------------------------------------
                0015 !  ---- public interfaces ----
                0016
                0017    public  lscale_cond, lscale_cond_init
                0018
                0019 !-----------------------------------------------------------------------
                0020 !   ---- version number ----
                0021
15388b052e Jean*0022  character(len=128) :: version = '$Id: lscale_cond_mod.F90,v 1.2 2017/08/11 20:48:51 jmc Exp $'
b2ea1d2979 Jean*0023  character(len=128) :: tag = '$Name:  $'
                0024
                0025 !-----------------------------------------------------------------------
                0026 !   ---- local/private data ----
                0027
                0028     real, parameter :: d622 = rdgas/rvgas
                0029     real, parameter :: d378 = 1.-d622
                0030
                0031     logical :: do_init=.true.
                0032
                0033 !-----------------------------------------------------------------------
                0034 !   --- namelist ----
                0035
                0036 real    :: hc=1.00
                0037 logical :: do_evap=.false.
                0038
                0039 namelist /lscale_cond_nml/  hc, do_evap
                0040
                0041 !-----------------------------------------------------------------------
                0042 !           description of namelist variables
                0043 !
                0044 !  hc        =  relative humidity at which large scale condensation
                0045 !               occurs, where 0 <= hc <= 1 (default: hc=1.)
                0046 !
                0047 !  do_evap   =  flag for the re-evaporation of moisture in
                0048 !               sub-saturated layers below, if do_evap=.true. then
                0049 !               re-evaporation is performed (default: do_evap=.false.)
                0050 !
                0051 !-----------------------------------------------------------------------
                0052
                0053 contains
                0054
                0055 !#######################################################################
                0056
                0057 !   subroutine lscale_cond (tin, qin, pfull, phalf, coldT, &
                0058 !                           rain, snow, tdel, qdel, mask, conv)
                0059 subroutine lscale_cond (tin, qin, pfull, phalf, coldT, &
                0060                         rain, snow, tdel, qdel, qsat,  &
                0061                         myThid, mask, conv )
                0062
                0063 !-----------------------------------------------------------------------
                0064 !
                0065 !                      large scale condensation
                0066 !
                0067 !-----------------------------------------------------------------------
                0068 !
                0069 !   input:  tin      temperature at full model levels
                0070 !           qin      specific humidity of water vapor at full
                0071 !                      model levels
                0072 !           pfull    pressure at full model levels
                0073 !           phalf    pressure at half (interface) model levels
                0074 !           coldT    should precipitation be snow at this point?
                0075 !   optional:
                0076 !           mask     optional mask (0 or 1.)
                0077 !           conv     logical flag; if true then no large-scale
                0078 !                       adjustment is performed at that grid-point or
                0079 !                       model level
                0080 !
                0081 !  output:  rain     liquid precipitation (kg/m2)
                0082 !           snow     frozen precipitation (kg/m2)
                0083 !           tdel     temperature tendency at full model levels
                0084 !           qdel     specific humidity tendency (of water vapor) at
                0085 !                      full model levels
                0086 !           qsat     saturated specific humidity
                0087 !
                0088 !-----------------------------------------------------------------------
                0089 !--------------------- interface arguments -----------------------------
                0090
                0091    real   , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf
                0092    logical   , intent(in) , dimension(:,:):: coldT
                0093    real   , intent(out), dimension(:,:)   :: rain,snow
                0094    real   , intent(out), dimension(:,:,:) :: tdel, qdel, qsat
                0095    integer, intent(in)                    :: myThid
                0096    real   , intent(in) , dimension(:,:,:), optional :: mask
                0097    logical, intent(in) , dimension(:,:,:), optional :: conv
                0098 !-----------------------------------------------------------------------
                0099 !---------------------- local data -------------------------------------
                0100
                0101 logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust
                0102    real,dimension(size(tin,1),size(tin,2),size(tin,3)) ::  &
                0103                              esat, desat, dqsat, pmes, pmass
                0104    real,dimension(size(tin,1),size(tin,2))             :: hlcp, precip
                0105 integer  k, kx
                0106 !-----------------------------------------------------------------------
                0107 !     computation of precipitation by condensation processes
                0108 !-----------------------------------------------------------------------
                0109
                0110 !      if (do_init) call error_mesg ('lscale_cond',  &
                0111 !                         'lscale_cond_init has not been called.', FATAL)
                0112
                0113       kx=size(tin,3)
                0114
                0115 !----- compute proper latent heat --------------------------------------
                0116       WHERE (coldT)
                0117            hlcp = HLs/Cp_air
                0118       ELSEWHERE
                0119            hlcp = HLv/Cp_air
                0120       END WHERE
                0121
                0122 !----- saturation vapor pressure (esat) & specific humidity (qsat) -----
                0123
                0124       call  escomp (tin,esat)
                0125       call descomp (tin,desat)
                0126
                0127       esat(:,:,:)=esat(:,:,:)*hc
                0128
                0129    where (pfull(:,:,:) > d378*esat(:,:,:))
                0130       pmes(:,:,:)=1.0/(pfull(:,:,:)-d378*esat(:,:,:))
                0131       qsat(:,:,:)=d622*esat(:,:,:)*pmes(:,:,:)
                0132       qsat(:,:,:)=max(0.0,qsat(:,:,:))
                0133      dqsat(:,:,:)=d622*pfull(:,:,:)*desat(:,:,:)*pmes(:,:,:)*pmes(:,:,:)
                0134    elsewhere
                0135       pmes(:,:,:)=0.0
                0136       qsat(:,:,:)=0.0
                0137      dqsat(:,:,:)=0.0
                0138    endwhere
                0139
                0140 !--------- do adjustment where greater than saturated value ------------
                0141
                0142    if (present(conv)) then
                0143 !!!!  do_adjust(:,:,:)=(.not.conv(:,:,:) .and. qin(:,:,:) > qsat(:,:,:))
                0144       do_adjust(:,:,:)=(.not.conv(:,:,:) .and.   &
                0145                          (qin(:,:,:) - qsat(:,:,:))*qsat(:,:,:) > 0.0)
                0146    else
                0147 !!!!  do_adjust(:,:,:)=(qin(:,:,:) > qsat(:,:,:))
                0148       do_adjust(:,:,:)=( (qin(:,:,:) - qsat(:,:,:))*qsat(:,:,:) > 0.0)
                0149    endif
                0150
                0151    if (present(mask)) then
                0152       do_adjust(:,:,:)=do_adjust(:,:,:) .and. (mask(:,:,:) > 0.5)
                0153    end if
                0154
                0155 !----------- compute adjustments to temp and spec humidity -------------
                0156    do k = 1,kx
                0157    where (do_adjust(:,:,k))
                0158       qdel(:,:,k)=(qsat(:,:,k)-qin(:,:,k))/(1.0+hlcp(:,:)*dqsat(:,:,k))
                0159       tdel(:,:,k)=-hlcp(:,:)*qdel(:,:,k)
                0160    elsewhere
                0161       qdel(:,:,k)=0.0
                0162       tdel(:,:,k)=0.0
                0163    endwhere
                0164    end do
                0165
                0166 !------------ pressure mass of each layer ------------------------------
                0167
                0168    do k=1,kx
                0169       pmass(:,:,k)=(phalf(:,:,k+1)-phalf(:,:,k))/Grav
                0170    enddo
                0171
                0172 !------------ re-evaporation of precipitation in dry layer below -------
                0173
                0174    if (do_evap) then
                0175       if (present(mask)) then
                0176 !        call precip_evap (pmass,tin,qin,qsat,dqsat,hlcp,tdel,qdel,mask)
                0177          call precip_evap (pmass,tin,qin,qsat,dqsat,hlcp,tdel,qdel,      &
                0178                            myThid, mask )
                0179       else
                0180          call precip_evap (pmass,tin,qin,qsat,dqsat,hlcp,tdel,qdel,      &
                0181                            myThid )
                0182       endif
                0183    endif
                0184
                0185 !------------ integrate precip -----------------------------------------
                0186
                0187       precip(:,:)=0.0
                0188    do k=1,kx
                0189       precip(:,:)=precip(:,:)-pmass(:,:,k)*qdel(:,:,k)
                0190    enddo
                0191       precip(:,:)=max(precip(:,:),0.0)
                0192
                0193    !assign precip to snow or rain
                0194    WHERE (coldT)
                0195       snow = precip
                0196       rain = 0.
                0197    ELSEWHERE
                0198       rain = precip
                0199       snow = 0.
                0200    END WHERE
                0201
                0202 !-----------------------------------------------------------------------
                0203
                0204    end subroutine lscale_cond
                0205
                0206 !#######################################################################
                0207
                0208 ! subroutine precip_evap (pmass, tin, qin, qsat, dqsat, hlcp, &
                0209 !                         tdel, qdel, mask)
                0210 subroutine precip_evap (pmass, tin, qin, qsat, dqsat, hlcp, &
                0211                         tdel, qdel, myThid, mask )
                0212
                0213 !-----------------------------------------------------------------------
                0214 !        performs re-evaporation of falling precipitation
                0215 !-----------------------------------------------------------------------
                0216    real, intent(in),    dimension(:,:,:) :: pmass, tin, qin, qsat, dqsat
                0217    real, intent(in),    dimension(:,:)   :: hlcp
                0218    real, intent(inout), dimension(:,:,:) :: tdel, qdel
                0219    integer, intent(in)                   :: myThid
                0220    real, intent(in), dimension(:,:,:), optional :: mask
                0221 !-----------------------------------------------------------------------
                0222    real, dimension(size(tin,1),size(tin,2)) :: exq, def
                0223
                0224    integer  k
                0225 !-----------------------------------------------------------------------
                0226     exq(:,:)=0.0
                0227
                0228     do k=1,size(tin,3)
                0229
                0230         where (qdel(:,:,k) < 0.0)  exq(:,:) = exq(:,:) -  &
                0231                                                qdel(:,:,k)*pmass(:,:,k)
                0232
                0233         if (present(mask)) exq(:,:) = exq(:,:)*mask(:,:,k)
                0234
                0235 !  ---- evaporate precip where needed ------
                0236
                0237         where ( (qdel(:,:,k) >= 0.0) .and. (exq(:,:) > 0.0) )
                0238             exq(:,:) = exq(:,:) / pmass(:,:,k)
                0239             def(:,:) = (qsat(:,:,k)-qin(:,:,k))/(1.+hlcp(:,:)*dqsat(:,:,k))
                0240             def(:,:) = min(max(def(:,:),0.0),exq(:,:))
                0241             qdel(:,:,k) = qdel(:,:,k) + def(:,:)
                0242             tdel(:,:,k) = tdel(:,:,k) - def(:,:)*hlcp(:,:)
                0243             exq(:,:) = (exq(:,:)-def(:,:))*pmass(:,:,k)
                0244         endwhere
                0245
                0246     enddo
                0247
                0248 !-----------------------------------------------------------------------
                0249
                0250    end subroutine precip_evap
                0251
                0252 !#######################################################################
                0253
                0254    subroutine lscale_cond_init (myThid)
                0255
                0256 !-----------------------------------------------------------------------
                0257 !
                0258 !        initialization for large scale condensation
                0259 !
                0260 !-----------------------------------------------------------------------
                0261
                0262   integer, intent(in) :: myThid
                0263
                0264 !-----------------------------------------------------------------------
                0265
                0266 ! integer  unit,io,ierr
                0267   integer         :: iUnit
                0268   CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0269
                0270 !----------- read namelist ---------------------------------------------
                0271
                0272 !    _BARRIER
                0273 !    _BEGIN_MASTER(myThid)
                0274      CALL BARRIER(myThid)
                0275      IF ( myThid.EQ.1 ) THEN
                0276
                0277      WRITE(msgBuf,'(A)') 'LSCALE_COND_INIT: opening data.atm_gray'
                0278      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0279      CALL OPEN_COPY_DATA_FILE(                                      &
                0280                            'data.atm_gray', 'LSCALE_COND_INIT',     &
                0281                            iUnit,                                   &
                0282                            myThid )
                0283 !    Read parameters from open data file
                0284      READ(UNIT=iUnit,NML=lscale_cond_nml)
                0285      WRITE(msgBuf,'(A)')                                            &
                0286           'LSCALE_COND_INIT: finished reading data.atm_gray'
                0287      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0288 !    Close the open data file
15388b052e Jean*0289 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0290      CLOSE(iUnit)
15388b052e Jean*0291 #else
                0292      CLOSE(iUnit,STATUS='DELETE')
                0293 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0294
                0295 !      if (file_exist('input.nml')) then
                0296 !         unit = open_file (file='input.nml', action='read')
                0297 !         ierr=1; do while (ierr /= 0)
                0298 !            read  (unit, nml=lscale_cond_nml, iostat=io, end=10)
                0299 !            ierr = check_nml_error (io,'lscale_cond_nml')
                0300 !         enddo
                0301 !  10     call close_file (unit)
                0302 !      endif
                0303
                0304 !---------- output namelist --------------------------------------------
                0305
                0306 !      unit = open_file (file='logfile.out', action='append')
                0307 !      if ( mpp_pe() == 0 ) then
                0308 !           write (unit,'(/,80("="),/(a))') trim(version), trim(tag)
                0309 !           write (unit,nml=lscale_cond_nml)
                0310 !      endif
                0311 !      call close_file (unit)
                0312
                0313       do_init=.false.
                0314
                0315      ENDIF
                0316      CALL BARRIER(myThid)
                0317
                0318    end subroutine lscale_cond_init
                0319
                0320 !#######################################################################
                0321
                0322 end module lscale_cond_mod
                0323