Newer
Older
!MNH_LIC Copyright 2006-2022 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
SUBROUTINE TH_R_FROM_THL_RT(CST, NEBN, KT, HFRAC_ICE,PFRAC_ICE,PP, &
PTHL, PRT, PTH, PRV, PRL, PRI, &
PRSATW, PRSATI, PRR, PRS, PRG, PRH, OOCEAN,&
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ******
! => Don't use drHook !!!
! "compute_frac_ice.func.h" must be included at the same time
! #################################################################
!
!
!!**** *TH_R_FROM_THL_RT* - computes the non-conservative variables
!! from conservative variables
!!
!! PURPOSE
!! -------
!!
!!** METHOD
!! ------
!!
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!!
!! REFERENCE
!! ---------
!!
!! AUTHOR
!! ------
!! Julien PERGAUD * Meteo-France *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/03/06
!! S. Riette April 2011 : ice added, allow ZRLTEMP to be negative
!! we use dQsat/dT to help convergence
!! use of optional PRR, PRS, PRG, PRH
!! S. Riette Nov 2016: support for HFRAC_ICE='S'
!!
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST, ONLY : CST_t
!
IMPLICIT NONE
!
!
!* 0.1 declarations of arguments
!
TYPE(CST_t), INTENT(IN) :: CST
INTEGER, INTENT(IN) :: KT
CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
LOGICAL, INTENT(IN) :: OOCEAN ! switch for Ocean model version
REAL, DIMENSION(KT), INTENT(INOUT) :: PFRAC_ICE
REAL, DIMENSION(KT), INTENT(IN) :: PP ! Pressure
REAL, DIMENSION(KT), INTENT(IN) :: PTHL ! thetal to transform into th
REAL, DIMENSION(KT),INTENT(IN) :: PRT ! Total mixing ratios to transform into rv,rc and ri
REAL, DIMENSION(KT),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
REAL, DIMENSION(KT), INTENT(OUT):: PTH ! th
REAL, DIMENSION(KT), INTENT(OUT):: PRV ! vapor mixing ratio
REAL, DIMENSION(KT), INTENT(INOUT):: PRL ! vapor mixing ratio
REAL, DIMENSION(KT), INTENT(INOUT):: PRI ! vapor mixing ratio
REAL, DIMENSION(KT), INTENT(OUT) :: PRSATW ! estimated mixing ration at saturation over water
REAL, DIMENSION(KT), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturation over ice
REAL, DIMENSION(KT, 16), INTENT(OUT) :: PBUF ! buffer to replace automatic arrays
INTEGER, OPTIONAL :: KB !first index to deal with (default is 1)
INTEGER, OPTIONAL :: KE !last index to deal with (default if KT)
!
!-------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
INTEGER :: II ! Loop control
INTEGER :: JITER ! number of iterations
INTEGER, PARAMETER :: IEXN=1, IRVSAT=2, ICPH=3, IRLTEMP=4, ICPH2=5, IT=6, ILVOCPEXN=7, ILSOCPEXN=8, &
& IDRSATODT=9, IDRSATODTW=10, IDRSATODTI=11, IFOESW=12, IFOESI=13, &
& ILOGT=14, I99PP=15, I1PRT=16
REAL :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
!----------------------------------------------------------------------------
!
!* 1 Initialisation
! --------------
!
!
!
IB=MERGE(KB, 1, PRESENT(KB))
IE=MERGE(KE, KT, PRESENT(KE))
!Number of iterations
JITER=2
!
!Computation of PBUF(IB:IE, ICPH2) depending on dummy arguments received
PBUF(IB:IE, ICPH2)=0
IF(PRESENT(PRR)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCL*PRR(IB:IE)
IF(PRESENT(PRS)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRS(IB:IE)
IF(PRESENT(PRG)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRG(IB:IE)
IF(PRESENT(PRH)) PBUF(IB:IE, ICPH2)=PBUF(IB:IE, ICPH2) + CST%XCI*PRH(IB:IE)
!
!Computation of an approximate state thanks to PRL and PRI guess
PBUF(IB:IE, IEXN)=(PP(IB:IE)/CST%XP00) ** CST%RDSCPD
PBUF(J, I99PP)=0.99*PP(J)
PRV(J)=PRT(J)-PRL(J)-PRI(J)
PBUF(J, ICPH)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + PBUF(J, ICPH2)
ZVAR2=PBUF(J, ICPH)*PBUF(J, IEXN)
ZDELT=(PTHL(J)*PBUF(J, IEXN))-CST%XTT
PBUF(J, ILVOCPEXN) = (CST%XLVTT + (CST%XCPV-CST%XCL) * ZDELT) /ZVAR2
PBUF(J, ILSOCPEXN) = (CST%XLSTT + (CST%XCPV-CST%XCI) * ZDELT) /ZVAR2
PTH(J)=PTHL(J)+PBUF(J, ILVOCPEXN)*PRL(J)+PBUF(J, ILSOCPEXN)*PRI(J)
PBUF(J, I1PRT)=1+PRT(J)
ENDDO
!
!
! 2 Iteration
! ---------
DO II=1,JITER
IF (OOCEAN) THEN
PBUF(IB:IE, IT)=PTH(IB:IE)
PBUF(IB:IE, IT)=PTH(IB:IE)*PBUF(IB:IE, IEXN)
END IF
!Computation of liquid/ice fractions
PFRAC_ICE(IB:IE) = 0.
DO J=IB, IE
IF(PRL(J)+PRI(J) > 1.E-20) THEN
PFRAC_ICE(J) = PRI(J) / (PRL(J)+PRI(J))
ENDIF
ENDDO
CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEBN,PFRAC_ICE(IB:IE),PBUF(IB:IE, IT))
!Computation of Rvsat and dRsat/dT
!In this version QSAT, QSATI, DQSAT and DQASATI functions are not used
!due to performance issue
! Log does not vectorize on all compilers:
PBUF(IB:IE, ILOGT)=LOG(PBUF(IB:IE, IT))
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
PBUF(J, IFOESW) = MIN(EXP( CST%XALPW - CST%XBETAW/PBUF(J, IT) - CST%XGAMW*PBUF(J, ILOGT) ), PBUF(J, I99PP))
PBUF(J, IFOESI) = MIN(EXP( CST%XALPI - CST%XBETAI/PBUF(J, IT) - CST%XGAMI*PBUF(J, ILOGT) ), PBUF(J, I99PP))
PRSATW(J) = CST%XRD/CST%XRV*PBUF(J, IFOESW)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESW)/PP(J))
PRSATI(J) = CST%XRD/CST%XRV*PBUF(J, IFOESI)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESI)/PP(J))
ZTPOW2=PBUF(J, IT)**2
PBUF(J, IDRSATODTW) = PRSATW(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESW)/PP(J) ) &
* (CST%XBETAW/ZTPOW2 - CST%XGAMW/PBUF(J, IT))*PBUF(J, I1PRT)
PBUF(J, IDRSATODTI) = PRSATI(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESI)/PP(J) ) &
* (CST%XBETAI/ZTPOW2 - CST%XGAMI/PBUF(J, IT))*PBUF(J, I1PRT)
!PRSATW(J) = QSAT(PBUF(J, IT),PP(J)) !qsatw
!PRSATI(J) = QSATI(PBUF(J, IT),PP(J)) !qsati
!PBUF(J, IDRSATODTW) = DQSAT(PBUF(J, IT),PP(J),PRSATW(J))*PBUF(J, I1PRT)
!PBUF(J, IDRSATODTI) = DQSATI(PBUF(J, IT),PP(J),PRSATI(J))*PBUF(J, I1PRT)
PRSATW(J) = PRSATW(J)*PBUF(J, I1PRT)
PRSATI(J) = PRSATI(J)*PBUF(J, I1PRT)
PBUF(J, IRVSAT) = PRSATW(J)*(1-PFRAC_ICE(J)) + PRSATI(J)*PFRAC_ICE(J)
PBUF(J, IDRSATODT) = (PBUF(J, IDRSATODTW)*(1-PFRAC_ICE(J))+ &
& PBUF(J, IDRSATODTI)*PFRAC_ICE(J))
!Computation of new PRL, PRI and PRV
!Correction term applied to (PRV(J)-PBUF(J, IRVSAT)) is computed assuming that
!PBUF(J, ILVOCPEXN), PBUF(J, ILSOCPEXN) and PBUF(J, ICPH) don't vary to much with T. It takes into account
!the variation (estimated linear) of Qsat with T
PBUF(J, IRLTEMP)=(PRV(J)-PBUF(J, IRVSAT))/ &
&(1 + PBUF(J, IDRSATODT)*PBUF(J, IEXN)* &
& (PBUF(J, ILVOCPEXN)*(1-PFRAC_ICE(J))+PBUF(J, ILSOCPEXN)*PFRAC_ICE(J)))
PBUF(J, IRLTEMP)=MIN(MAX(-PRL(J)-PRI(J), PBUF(J, IRLTEMP)),PRV(J))
PRV(J)=PRV(J)-PBUF(J, IRLTEMP)
PRL(J)=PRL(J)+PRI(J)+PBUF(J, IRLTEMP)
PRI(J)=PFRAC_ICE(J) * (PRL(J))
PRL(J)=(1-PFRAC_ICE(J)) * (PRT(J) - PRV(J))
!Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
PBUF(J, ICPH)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + PBUF(J, ICPH2)
!Computation of L/Cph/EXN, then new PTH
ZVAR2=PBUF(J, ICPH)*PBUF(J, IEXN)
PBUF(J, ILVOCPEXN) = (CST%XLVTT + (CST%XCPV-CST%XCL) * (PBUF(J, IT)-CST%XTT)) /ZVAR2
PBUF(J, ILSOCPEXN) = (CST%XLSTT + (CST%XCPV-CST%XCI) * (PBUF(J, IT)-CST%XTT)) /ZVAR2
PTH(J)=PTHL(J)+PBUF(J, ILVOCPEXN)*PRL(J)+PBUF(J, ILSOCPEXN)*PRI(J)
!Computation of estimated mixing ration at saturation
!To compute the adjustement a first order development was used
ZVAR1=PTH(J)*PBUF(J, IEXN)-PBUF(J, IT)
PRSATW(J)=PRSATW(J) + PBUF(J, IDRSATODTW)*ZVAR1
PRSATI(J)=PRSATI(J) + PBUF(J, IDRSATODTI)*ZVAR1
ENDDO
ENDDO
END SUBROUTINE TH_R_FROM_THL_RT