!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module moist_conv_mod !----------------------------------------------------------------------- use mpp_mod, only : mpp_pe, & mpp_root_pe, & stdlog use time_manager_mod, only : time_type use Diag_Manager_Mod, ONLY: register_diag_field, send_data use sat_vapor_pres_mod, ONLY: lookup_es_des, compute_qs, descomp use fms_mod, ONLY: error_mesg, file_exist, open_namelist_file, & check_nml_error, close_file, & FATAL, WARNING, NOTE, mpp_pe, mpp_root_pe, & write_version_number, stdlog use constants_mod, ONLY: HLv, HLs, cp_air, grav, rdgas, rvgas use fms_mod, only : write_version_number, ERROR_MESG, FATAL use field_manager_mod, only : MODEL_ATMOS use tracer_manager_mod, only : get_tracer_index, & get_number_tracers, & get_tracer_names, & get_tracer_indices, & query_method, & NO_TRACER implicit none private !------- interfaces in this module ------------ public :: moist_conv, moist_conv_Init, moist_conv_end !----------------------------------------------------------------------- !---- namelist ---- real :: HC = 1.00 real :: beta = 0.0 real :: TOLmin=.02, TOLmax=.10 integer :: ITSMOD=30 logical :: do_simple =.false. !----- note beta is the fraction of convective condensation that is !----- detrained into a stratiform cloud namelist /moist_conv_nml/ HC, beta, TOLmin, TOLmax, ITSMOD, do_simple !----------------------------------------------------------------------- !---- VERSION NUMBER ----- character(len=128) :: version = '$Id: moist_conv.F90,v 17.0.4.1 2009/09/24 18:23:42 wfc Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: module_is_initialized = .false. !---------- initialize constants used by this module ------------------- real, parameter :: d622 = rdgas/rvgas real, parameter :: d378 = 1.0-d622 real, parameter :: grav_inv = 1.0/grav real, parameter :: rocp = rdgas/cp_air real :: missing_value = -999. integer :: nsphum, nql, nqi, nqa ! tracer indices for stratiform clouds integer :: id_tdt_conv, id_qdt_conv, id_prec_conv, id_snow_conv, & id_qldt_conv, id_qidt_conv, id_qadt_conv, & id_ql_conv_col, id_qi_conv_col, id_qa_conv_col,& id_q_conv_col, id_t_conv_col character(len=3) :: mod_name = 'mca' logical :: do_mca_tracer = .false. integer :: num_mca_tracers = 0 integer :: num_tracers integer, allocatable, dimension(:) :: id_tracer_conv, id_tracer_conv_col !----------------------------------------------------------------------- !----------------------------------------------------------------------- CONTAINS !####################################################################### subroutine moist_conv ( Tin, Qin, Pfull, Phalf, coldT, & ! required Tdel, Qdel, Rain, Snow, & ! required dtinv, Time, is, js, tracers, qtrmca, & ! required Lbot, mask, Conv, & ! optional ql, qi, cf, qldel, qidel, cfdel) ! optional !----------------------------------------------------------------------- ! ! MOIST CONVECTIVE ADJUSTMENT ! !----------------------------------------------------------------------- ! ! INPUT: Tin temperature at full model levels ! Qin specific humidity of water vapor at full ! model levels ! Pfull pressure at full model levels ! Phalf pressure at half model levels ! coldT Should MCA produce snow in this column? ! ! OUTPUT: Tdel temperature adjustment at full model levels (deg k) ! Qdel specific humidity adjustment of water vapor at ! full model levels ! Rain liquid precipitiation (in Kg m-2) ! Snow ice phase precipitation (kg m-2) ! OPTIONAL ! ! INPUT: Lbot integer index of the lowest model level, ! Lbot is always <= size(Tin,3) ! ql liquid water condensate ! qi ice condensate ! cf stratiform cloud fraction (used only when ! operating with stratiform cloud scheme) (fraction) ! ! OUTPUT: Conv logical flag; TRUE then moist convective ! adjustment was performed at that model level. ! cfdel change in stratiform cloud fraction (fraction) ! qldel change in liquid water condensate due to ! convective detrainment (kg condensate /kg air) ! qidel change in ice condensate due to ! convective detrainment (kg condensate /kg air) ! !----------------------------------------------------------------------- !----------------------PUBLIC INTERFACE ARRAYS-------------------------- real, intent(INOUT), dimension(:,:,:) :: Tin, Qin real, intent(IN) , dimension(:,:,:) :: Pfull, Phalf logical, intent(IN) , dimension(:,:) :: coldT real, intent(OUT), dimension(:,:,:) :: Tdel, Qdel real, intent(OUT), dimension(:,:) :: Rain, Snow real, intent(IN) :: dtinv type(time_type), intent(in) :: Time integer, intent(IN) :: is, js real, dimension(:,:,:,:), intent(in) :: tracers real, dimension(:,:,:,:), intent(out) :: qtrmca integer, intent(IN) , dimension(:,:), optional :: Lbot real, intent(IN) , dimension(:,:,:), optional :: mask logical, intent(OUT), dimension(:,:,:), optional :: Conv real, intent(INOUT), dimension(:,:,:), optional :: ql, qi, cf real, intent(OUT), dimension(:,:,:), optional :: qldel, qidel, cfdel !----------------------------------------------------------------------- !----------------------PRIVATE (LOCAL) ARRAYS--------------------------- ! logical, dimension(size(Tin,1),size(Tin,2),size(Tin,3)) :: DO_ADJUST ! real, dimension(size(Tin,1),size(Tin,2),size(Tin,3)) :: & !----------------------------------------------------------------------- !----------------------PRIVATE (LOCAL) ARRAYS--------------------------- integer, dimension(size(Tin,1),size(Tin,2)) :: ISMVF integer, dimension(size(Tin,1),size(Tin,2),size(Tin,3)) :: IVF real, dimension(size(Tin,1),size(Tin,2),size(Tin,3)) :: & Qdif,Temp,Qmix,Esat,Qsat,Test1,Test2, & Esdiff_v real, dimension(size(Tin,1),size(Tin,2),size(Tin,3)-1) :: & Thalf,DelPoP,Esm,Esd,ALRM real, dimension(size(Tin,3)) :: C,Ta,Qa real, dimension(size(Tin,1),size(Tin,2)) :: HL integer :: i,j,k,kk,KX,ITER,MXLEV,MXLEV1,kstart,KTOP,KBOT,KBOTM1 real :: ALTOL,Sum0,Sum1,Sum2,EsDiff,EsVal,Thaf,Pdelta logical :: cloud_tracers_present real, dimension(size(Phalf,1),size(Phalf,2),size(Phalf,3)) :: pmass real, dimension(size(Phalf,1),size(Phalf,2)) :: tempdiag integer :: tr, num_cld_tracers logical :: used !----------------------------------------------------------------------- if (.not. module_is_initialized) call ERROR_MESG( 'MCA', & 'moist_conv_init has not been called', FATAL ) num_cld_tracers = count( (/present(ql),present(qi),present(cf),present(qldel),present(qidel),present(cfdel)/) ) if(num_cld_tracers == 0) then cloud_tracers_present = .false. else if(num_cld_tracers == 6) then cloud_tracers_present = .true. else call error_mesg('moist_conv','Either all or none of the cloud tracers and their tendencies must be present',FATAL) endif do k=1,size(Tin,3) pmass(:,:,k) = (Phalf(:,:,k+1)-Phalf(:,:,k))/GRAV end do KX=size(Tin,3) !------ compute Proper HL if(do_simple) then HL = HLv else WHERE (coldT) HL = HLs ELSEWHERE HL = HLv END WHERE endif !------ convert spec hum to mixing ratio ------ Temp(:,:,:)=Tin(:,:,:) Qmix(:,:,:)=Qin(:,:,:) do k=1,KX-1 DelPoP(:,:,k)=(Pfull(:,:,k+1)-Pfull(:,:,k))/Phalf(:,:,k+1) enddo !-------------SATURATION VAPOR PRESSURE FROM ETABL---------------------- ! compute qs; also return dqsdT call compute_qs (Temp, Pfull, Qsat, dqsdT=Esdiff_v, hc = hc, & esat=Esat) Qdif(:,:,:)=Max(0.0,Qmix(:,:,:)-Qsat(:,:,:)) !----------------------------------------------------------------------- ! MOIST CONVECTIVE ADJUSTMENT !----------------------------------------------------------------------- ! *** Set initial tolerance *** ALTOL = TOLmin do k=1,KX-1 Thalf(:,:,k)=0.50*(Temp(:,:,k)+Temp(:,:,k+1)) enddo call lookup_es_des (Thalf, Esm, Esd) do k=1,KX-1 ALRM(:,:,k)=rocp*DelPoP(:,:,k)*Thalf(:,:,k) & *(Phalf(:,:,k+1)+d622*HL(:,:)*Esm(:,:,k)/Thalf(:,:,k)/rdgas) & /(Phalf(:,:,k+1)+d622*HL(:,:)*Esd(:,:,k)/cp_air) enddo IVF (:,:,KX)=0 Test1(:,:,KX)=0.0 Test2(:,:,KX)=0.0 do k=1,KX-1 Test1(:,:,k)=Temp(:,:,k+1)-Temp(:,:,k) Test2(:,:,k)=ALRM(:,:,k)+ALTOL-Test1(:,:,k) enddo !!!!! Test1(:,:,:)=0.0-Qdif(:,:,:) Test1(:,:,:)=(0.0-Qdif(:,:,:))*Qsat(:,:,:) !-------IVF=1 in unstable layers where both levels are saturated-------- do k=1,KX-1 where (Test1(:,:,k) < 0.0 .and. Test1(:,:,k+1) < 0.0 .and. & Test2(:,:,k) < 0.0) IVF(:,:,k)=1 elsewhere IVF(:,:,k)=0 endwhere enddo ! ------ Set convection flag (for optional output only) -------- if (Present(Conv)) then Conv(:,:,1)=(IVF(:,:,1) == 1) do k=1,KX-1 Conv(:,:,k+1)=(IVF(:,:,k) == 1 .or. IVF(:,:,k+1) == 1) enddo endif ! ----- Set counter for each column ----- ISMVF(:,:)=0 do k=1,KX-1 ISMVF(:,:)=ISMVF(:,:)+IVF(:,:,k) enddo !----------------------------------------------------------------------- !---------------LOOP OVER EACH VERTICAL COLUMN-------------------------- do j=1,size(Tin,2) OUTER_LOOP: do i=1,size(Tin,1) !----------------------------------------------------------------------- if (ISMVF(i,j) == 0) CYCLE if (Present(Lbot)) then MXLEV=Lbot(i,j) else MXLEV=KX endif MXLEV1=MXLEV-1 ! *** Re-set initial tolerance *** ALTOL = TOLmin !----------(return here after increasing tolerance)-------------------- 1450 CONTINUE !--------------Iterations at the same tolerance------------------------- do 1740 ITER=1,ITSMOD !----------------------------------------------------------------------- kstart=1 1500 CONTINUE !-------------TEST TO DETERMINE UNSTABLE LAYER BLOCKS------------------- !-------Find top (KTOP) and bottom (KBOT) of unstable layers------------ do k=kstart,MXLEV1 if (IVF(i,j,k) == 1) GO TO 1505 enddo CYCLE OUTER_LOOP 1505 KTOP=k do k=KTOP,MXLEV1 if (IVF(i,j,k+1) == 0) then KBOT=k+1 GO TO 1510 endif enddo KBOT=MXLEV 1510 CONTINUE !----------------------------------------------------------------------- KBOTM1=KBOT - 1 Sum1=0.0 Sum2=0.0 !----------------------------------------------------------------------- do 1630 k=KTOP,KBOT !----------------------------------------------------------------------- if(do_simple) then call DEsComp (Temp(i,j,k),EsDiff) C(k)=d622*HC*EsDiff/Pfull(i,j,k) else C(k)=Pfull(i,j,k)-d378*Esat(i,j,k) if (C(k) <= 0.0) then C(k)=0.0 else C(k) = esdiff_v(i,j,k) endif endif Sum0=0.0 if (k == KBOT) GO TO 1625 kk=k 1620 if (kk > KBOTM1) GO TO 1625 Sum0=Sum0+ALRM(i,j,kk) kk=kk+1 GO TO 1620 1625 CONTINUE Pdelta=Phalf(i,j,k+1)-Phalf(i,j,k) Sum1=Sum1 + Pdelta*((cp_air+HL(i,j)*C(k))*(Temp(i,j,k)+Sum0)+ & HL(i,j)*(Qmix(i,j,k)-Qsat(i,j,k))) Sum2=Sum2 + Pdelta*(cp_air+HL(i,j)*C(k)) !----------------------------------------------------------------------- 1630 CONTINUE !----------------------------------------------------------------------- Ta(KBOT)=Sum1/Sum2 k=KTOP 1645 if (k > KBOTM1) GO TO 1641 Sum0=0.0 kk=k 1640 if (kk > KBOTM1) GO TO 1642 Sum0=Sum0+ALRM(i,j,kk) kk=kk+1 GO TO 1640 1642 Ta(k)=Ta(KBOT)-Sum0 k=k+1 GO TO 1645 !---------UPDATE T,R,ES,Esm,Esd & Qsat FOR THE ADJUSTED POINTS---------- 1641 do k=KTOP,KBOT Qa(k)=Qsat(i,j,k)+C(k)*(Ta(k)-Temp(i,j,k)) Temp(i,j,k)=Ta(k) Qmix(i,j,k)=Qa(k) call compute_qs ( Temp(i,j,k), Pfull(i,j,k), Qsat(i,j,k), & dqsdT = Esdiff_v(i,j,k), hc =hc, esat= Esat(i,j,k)) Qdif(i,j,k)=Max(0.0,Qmix(i,j,k)-Qsat(i,j,k)) enddo do k=KTOP,KBOTM1 Thaf=0.50*(Temp(i,j,k)+Temp(i,j,k+1)) !DIR$ INLINE call lookup_es_des (Thaf, EsVal, Esdiff) !DIR$ NOINLINE Esm (i,j,k)=HC*EsVal Esd (i,j,k)=HC*EsDiff ALRM(i,j,k)=rocp*DelPoP(i,j,k)*Thaf* & (Phalf(i,j,k+1)+d622*HL(i,j)*Esm(i,j,k)/Thaf/rdgas)/ & (Phalf(i,j,k+1)+d622*HL(i,j)*Esd(i,j,k)/cp_air) enddo !------------Is this the bottom of the current column ???--------------- kstart=KBOT+1 if (kstart <= MXLEV1) GO TO 1500 !----------------------------------------------------------------------- if (ITER == ITSMOD) GO TO 1740 !----------------------------------------------------------------------- do k=1,MXLEV1 Thaf=0.50*(Temp(i,j,k)+Temp(i,j,k+1)) !DIR$ INLINE call lookup_es_des (Thaf, EsVal, EsDiff) !DIR$ NOINLINE Esm (i,j,k)=HC*EsVal Esd (i,j,k)=HC*EsDiff ALRM(i,j,k)=rocp*DelPoP(i,j,k)*Thaf* & (Phalf(i,j,k+1)+d622*HL(i,j)*Esm(i,j,k)/Thaf/rdgas)/ & (Phalf(i,j,k+1)+d622*HL(i,j)*Esd(i,j,k)/cp_air) enddo do k=1,MXLEV1 IVF(i,j,k)=0 !!!! if (Qdif(i,j,k) > 0.0 .and. Qdif(i,j,k+1) > 0.0 .and. & if (Qdif(i,j,k)*Qsat(i,j,k) > 0.0 .and. & Qdif(i,j,k+1)*Qsat(i,j,k+1) > 0.0 .and. & (Temp(i,j,k+1)-Temp(i,j,k)) > (ALRM(i,j,k)+ALTOL)) then IVF(i,j,k) = 1 endif enddo ! ------ reset optional convection flag ------ if (Present(Conv)) then Conv(i,j,1)=(IVF(i,j,1) == 1) do k=1,MXLEV1 Conv(i,j,k+1)=(IVF(i,j,k) == 1 .or. IVF(i,j,k+1) == 1) enddo endif ! ------ Are all layers sufficiently stable ??? ------ ISMVF(i,j)=0 do k=1,MXLEV1 ISMVF(i,j)=ISMVF(i,j)+IVF(i,j,k) enddo if (ISMVF(i,j) == 0) CYCLE OUTER_LOOP !----------------------------------------------------------------------- 1740 CONTINUE !----------------------------------------------------------------------- !---------Maximum iterations reached: Increase tolerance (ALTOL)-------- ALTOL=2.0*ALTOL !del WRITE (*,9902) I,ALTOL call error_mesg ('moist_conv', 'Tolerence (ALTOL) doubled', NOTE) if (ALTOL <= TOLmax) GO TO 1450 ! WRITE (*,9903) ! WRITE (*,9904) (k,Temp(i,j,k),Qmix(i,j,k),Qsat(i,j,k), & ! Qdif(i,j,k),ALRM(i,j,k),k=1,MXLEV1) ! WRITE (*,9904) (k,Temp(i,j,k),Qmix(i,j,k),Qsat(i,j,k), & ! Qdif(i,j,k) ,k=MXLEV,MXLEV) call error_mesg ('moist_conv', 'maximum iterations reached', WARNING) !----------------------------------------------------------------------- enddo OUTER_LOOP enddo !----------------------------------------------------------------------- !---------------------- END OF i,j LOOP -------------------------------- !----------------------------------------------------------------------- !---- call Convective Detrainment subroutine ----- if (cloud_tracers_present) then !reset quantities cfdel = 0. qldel = 0. qidel = 0. CALL CONV_DETR(Qmix,Qin,Phalf,Temp,cf,coldT,cfdel,qldel,qidel) endif !----- compute adjustments to temp and spec hum ---- Tdel(:,:,:)=Temp(:,:,:)-Tin(:,:,:) Qdel(:,:,:)=Qmix(:,:,:)-Qin(:,:,:) !----- integrate precip ----- Rain(:,:)=0.0 Snow(:,:)=0.0 do k =1,KX if(do_simple) then Rain(:,:)=Rain(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & Qdel(:,:,k)*grav_inv if (cloud_tracers_present) then Rain(:,:)=Rain(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & qldel(:,:,k)*grav_inv endif else WHERE(coldT(:,:)) Snow(:,:)=Snow(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & Qdel(:,:,k)*grav_inv ELSEWHERE Rain(:,:)=Rain(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & Qdel(:,:,k)*grav_inv END WHERE !subtract off detrained condensate from surface precip if (cloud_tracers_present) then WHERE(coldT(:,:)) Snow(:,:)=Snow(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & qidel(:,:,k)*grav_inv ELSEWHERE Rain(:,:)=Rain(:,:)+(Phalf(:,:,k)-Phalf(:,:,k+1))* & qldel(:,:,k)*grav_inv END WHERE end if endif enddo Rain(:,:)=Max(Rain(:,:),0.0) Snow(:,:)=Max(Snow(:,:),0.0) !----------------------------------------------------------------------- !----------------- PRINT FORMATS ----------------------------------- 9902 FORMAT( ' *** ALTOL DOUBLED IN CONVAD AT I=', & I5,' ,ALTOL=', F10.4 ) 9903 FORMAT(/,' *** DIVERGENCE IN MOIST CONVECTIVE ADJUSTMENT ',/, & 4X,'K',14X,'T',14X,'R',13X,'Qsat',14X,'Qdif',12X,'ALRM',/) 9904 FORMAT (I5,5E15.7) !----------------------------------------------------------------------- !------- update input values and compute tendency ------- Tin=Tin+Tdel; Qin=Qin+Qdel Tdel=Tdel*dtinv; Qdel=Qdel*dtinv Rain=Rain*dtinv; Snow=Snow*dtinv !------- update input values , compute and add on tendency ----------- !------- in the case of strat ----------- if (cloud_tracers_present) then ql(:,:,:)=ql(:,:,:)+qldel(:,:,:) qi(:,:,:)=qi(:,:,:)+qidel(:,:,:) cf(:,:,:)=cf(:,:,:)+cfdel(:,:,:) qldel(:,:,:)=qldel(:,:,:)*dtinv qidel(:,:,:)=qidel(:,:,:)*dtinv cfdel(:,:,:)=cfdel(:,:,:)*dtinv endif !--------------------------------------------------------------------- ! define the effect of moist convective adjustment on the tracer ! fields. code to do so does not currently exist. !--------------------------------------------------------------------- qtrmca = 0. !------- diagnostics for dt/dt_ras ------- if ( id_tdt_conv > 0 ) then used = send_data ( id_tdt_conv, Tdel, Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for dq/dt_ras ------- if ( id_qdt_conv > 0 ) then used = send_data ( id_qdt_conv, Qdel, Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for precip_ras ------- if ( id_prec_conv > 0 ) then used = send_data ( id_prec_conv, Rain+Snow, Time, is, js ) endif !------- diagnostics for snow_ras ------- if ( id_snow_conv > 0 ) then used = send_data ( id_snow_conv, Snow, Time, is, js ) endif !------- diagnostics for water vapor path tendency ---------- if ( id_q_conv_col > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + Qdel(:,:,k)*pmass(:,:,k) end do used = send_data ( id_q_conv_col, tempdiag, Time, is, js ) end if !------- diagnostics for dry static energy tendency --------- if ( id_t_conv_col > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + Tdel(:,:,k)*cp_air*pmass(:,:,k) end do used = send_data ( id_t_conv_col, tempdiag, Time, is, js ) end if !------- stratiform cloud tendencies from cumulus convection ------------ if (cloud_tracers_present) then !------- diagnostics for dql/dt from RAS or donner ------- if ( id_qldt_conv > 0 ) then used = send_data ( id_qldt_conv, qldel(:,:,:), Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for dqi/dt from RAS or donner ------- if ( id_qidt_conv > 0 ) then used = send_data ( id_qidt_conv, qidel(:,:,:), Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for dqa/dt from RAS or donner ------- if ( id_qadt_conv > 0 ) then used = send_data ( id_qadt_conv, cfdel(:,:,:), Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for liquid water path tendency ------ if ( id_ql_conv_col > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + qldel(:,:,k)*pmass(:,:,k) end do used = send_data ( id_ql_conv_col, tempdiag, Time, is, js ) end if !------- diagnostics for ice water path tendency --------- if ( id_qi_conv_col > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + qidel(:,:,k)*pmass(:,:,k) end do used = send_data ( id_qi_conv_col, tempdiag, Time, is, js ) end if !---- diagnostics for column integrated cloud mass tendency --- if ( id_qa_conv_col > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + cfdel(:,:,k)*pmass(:,:,k) end do used = send_data ( id_qa_conv_col, tempdiag, Time, is, js ) end if end if ! if ( cloud_tracers_present ) do tr = 1, num_mca_tracers !------- diagnostics for dtracer/dt from RAS ------------- if ( id_tracer_conv(tr) > 0 ) then used = send_data ( id_tracer_conv(tr), qtrmca(:,:,:,tr), Time, is, js, 1, & rmask=mask ) endif !------- diagnostics for column tracer path tendency ----- if ( id_tracer_conv_col(tr) > 0 ) then tempdiag(:,:)=0. do k=1,kx tempdiag(:,:) = tempdiag(:,:) + qtrmca(:,:,k,tr)*pmass(:,:,k) end do used = send_data ( id_tracer_conv_col(tr), tempdiag, Time, is, js ) end if enddo end subroutine moist_conv !####################################################################### !####################################################################### SUBROUTINE CONV_DETR(qvout,qvin,phalf,T,cf,coldT,cfdel,qldel,qidel) IMPLICIT NONE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! This subroutine takes a fraction of the water condensed ! by the convection scheme and detrains in the top level ! undergoing convective adjustment. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! VARIABLES ! ! ------ ! INPUT: ! ------ ! ! qvout water vapor specific humidity AFTER adjustment ! (kg vapor/kg air) ! qvin water vapor specific humidity BEFORE adjustment ! (kg vapor/kg air) ! phalf pressure at model half levels (Pascals) ! T Temperature (Kelvin) ! cf cloud fraction (fraction) ! coldT is condensation of ice nature? ! ! ------------- ! INPUT/OUTPUT: ! ------------- ! ! cfdel Change in cloud fraction due to detrainment (fraction) ! qldel Increase in liquid water due to detrainment ! (kg condensate/kg air) ! qidel Increase in ice due to detrainment ! (kg condensate/kg air) ! ! ------------------- ! INTERNAL VARIABLES: ! ------------------- ! ! precipsource accumulated source of precipitation ! (kg condensate /meter/ (seconds*squared)) ! ktop integer of top level undergoing convection ! accum logical variable indicating whether or not to ! add precip ! fT fraction of condensate that is liquid ! i,j,k looping variables ! IDIM,JDIM,KDIM dimensions of input arrays ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! User Interface variables ! ------------------------ REAL, INTENT (IN), DIMENSION(:,:,:) :: qvout,qvin,T,cf REAL, INTENT (IN), DIMENSION(:,:,:) :: phalf LOGICAL, INTENT (IN), DIMENSION(:,:) :: coldT REAL, INTENT (INOUT),DIMENSION(:,:,:):: cfdel,qldel,qidel ! Internal variables ! ------------------ INTEGER :: i,j,k,IDIM,JDIM,KDIM,ktop LOGICAL :: accum REAL :: precipsource ! ! Code ! ---- ! reinitialize variables cfdel(:,:,:) = 0. qidel(:,:,:) = 0. qldel(:,:,:) = 0. IDIM = SIZE(qvout,1) JDIM = SIZE(qvout,2) KDIM = SIZE(qvout,3) !---loop over grid columns----! DO i = 1, IDIM DO j = 1, JDIM !reset variables precipsource = 0. accum = .FALSE. DO k = 1, KDIM !begin new convective event IF ((qvout(i,j,k) .ne. qvin(i,j,k)) .and. & (.NOT. accum)) THEN ktop = k accum = .TRUE. END IF !if convective event is over compute detrainment IF ( (accum) .and. (qvout(i,j,k) .eq. qvin(i,j,k)) & .and. (precipsource .gt. 0.)) THEN if (coldT(i,j)) then qidel(i,j,ktop) = beta * precipsource / & (phalf(i,j,ktop+1)-phalf(i,j,ktop)) else qldel(i,j,ktop) = beta * precipsource / & (phalf(i,j,ktop+1)-phalf(i,j,ktop)) end if cfdel(i,j,ktop) = MAX(0.,HC-cf(i,j,ktop)) accum = .FALSE. precipsource = 0. END IF !accumulate precip IF (accum) THEN precipsource = precipsource + & ( qvin(i,j,k) -qvout(i,j,k))* & (phalf(i,j,k+1)-phalf(i,j,k)) END IF END DO !---end k loop over vertical column !---clear any remaining precip IF ( (precipsource .gt. 0.) .and. (accum) ) THEN if (coldT(i,j)) then qidel(i,j,ktop) = beta * precipsource / & (phalf(i,j,ktop+1)-phalf(i,j,ktop)) else qldel(i,j,ktop) = beta * precipsource / & (phalf(i,j,ktop+1)-phalf(i,j,ktop)) end if END IF END DO !---end j loop END DO !---end i loop END SUBROUTINE CONV_DETR !####################################################################### subroutine moist_conv_init (axes, Time, tracers_in_mca) integer, intent(in) :: axes(4) type(time_type), intent(in) :: Time logical, dimension(:), intent(in), optional :: tracers_in_mca !----------------------------------------------------------------------- integer :: unit, io, ierr, logunit integer :: nn, tr character(len=128) :: diagname, diaglname, tendunits, name, units !----------------------------------------------------------------------- if (file_exist('input.nml')) then unit = open_namelist_file () ierr=1; do while (ierr /= 0) read (unit, nml=moist_conv_nml, iostat=io, end=10) ierr = check_nml_error (io,'moist_conv_nml') enddo 10 call close_file (unit) endif !---------- output namelist -------------------------------------------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number (version, tagname) logunit = stdlog() write (logunit,nml=moist_conv_nml) endif id_tdt_conv = register_diag_field ( mod_name, & 'tdt_conv', axes(1:3), Time, & 'Temperature tendency from moist conv adj', 'deg_K/s', & missing_value=missing_value ) id_qdt_conv = register_diag_field ( mod_name, & 'qdt_conv', axes(1:3), Time, & 'Spec humidity tendency from moist conv adj', 'kg/kg/s', & missing_value=missing_value ) id_prec_conv = register_diag_field ( mod_name, & 'prec_conv', axes(1:2), Time, & 'Precipitation rate from moist conv adj', 'kg/m2/s' ) id_snow_conv = register_diag_field ( mod_name, & 'snow_conv', axes(1:2), Time, & 'Frozen precip rate from moist conv adj', 'kg/m2/s' ) id_q_conv_col = register_diag_field ( mod_name, & 'q_conv_col', axes(1:2), Time, & 'Water vapor path tendency from moist conv adj','kg/m2/s' ) id_t_conv_col = register_diag_field ( mod_name, & 't_conv_col', axes(1:2), Time, & 'Column static energy tendency from moist conv adj','W/m2' ) !--------------------------------------------------------------------- ! --- Find the tracer indices !--------------------------------------------------------------------- call get_number_tracers (MODEL_ATMOS, num_prog= num_tracers) if ( num_tracers .gt. 0 ) then else call error_mesg('moist_conv_init', 'No atmospheric tracers found', FATAL) endif ! get tracer indices for stratiform cloud variables nsphum = get_tracer_index ( MODEL_ATMOS, 'sphum' ) nql = get_tracer_index ( MODEL_ATMOS, 'liq_wat' ) nqi = get_tracer_index ( MODEL_ATMOS, 'ice_wat' ) nqa = get_tracer_index ( MODEL_ATMOS, 'cld_amt' ) !---------------------------------------------------------------------- ! determine how many tracers are to be transported by moist_conv_mod. !---------------------------------------------------------------------- num_mca_tracers = count(tracers_in_mca) if (num_mca_tracers > 0) then do_mca_tracer = .true. else do_mca_tracer = .false. endif !--------------------------------------------------------------------- ! allocate the arrays to hold the diagnostics for the moist_conv ! tracers. !--------------------------------------------------------------------- allocate(id_tracer_conv (num_mca_tracers)) ; id_tracer_conv = 0 allocate(id_tracer_conv_col(num_mca_tracers)) ; id_tracer_conv_col = 0 nn = 1 do tr = 1,num_tracers if (tracers_in_mca(tr)) then call get_tracer_names(MODEL_ATMOS, tr, name=name, units=units) !---------------------------------------------------------------------- ! for the column tendencies, the name for the diagnostic will be ! the name of the tracer followed by 'dt_MCA'. the longname will be ! the name of the tracer followed by ' tendency from MCA'. units are ! the supplied units of the tracer divided by seconds. !---------------------------------------------------------------------- diagname = trim(name)//'dt_MCA' diaglname = trim(name)//' tendency from MCA' tendunits = trim(units)//'/s' id_tracer_conv(nn) = register_diag_field ( mod_name, & trim(diagname), axes(1:3), Time, & trim(diaglname), trim(tendunits), & missing_value=missing_value ) !---------------------------------------------------------------------- ! for the column integral tendencies, the name for the diagnostic ! will be the name of the tracer followed by 'dt_MCA_col'. the long- ! name will be the name of the tracer followed by ' path tendency ! from MCA'. units are the supplied units of the tracer multiplied ! by m**2 /kg divided by seconds. !---------------------------------------------------------------------- diagname = trim(name)//'dt_MCA_col' diaglname = trim(name)//' path tendency from MCA' tendunits = trim(units)//'m2/kg/s' id_tracer_conv_col(nn) = register_diag_field ( mod_name, & trim(diagname), axes(1:2), Time, & trim(diaglname), trim(tendunits), & missing_value=missing_value) nn = nn + 1 endif end do module_is_initialized = .true. !----------------------------------------------------------------------- end subroutine moist_conv_init !####################################################################### subroutine moist_conv_end integer :: log_unit if(.not.module_is_initialized) then return else module_is_initialized = .FALSE. endif log_unit = stdlog() if ( mpp_pe() == mpp_root_pe() ) then write (log_unit,'(/,(a))') 'Exiting moist_conv.' endif end subroutine moist_conv_end !####################################################################### end module moist_conv_mod