diff --git a/docs/TODO b/docs/TODO
index db758901a11f1e5837793b6ced157f19213930e1..00d21af05f14998f9fe9f1cdfdc1030fedc92379 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -68,8 +68,10 @@ Nettoyage apl_arome non fait (pb a la compilation) ==> 4 arguments dans aro_turb
turb.F90 : il reste un CALL à SOURCES_NEG_CORRECT à ajouter. Besoin de récupérer CCLOUD dans apl_arome : comment ?
-Regarder s'il ne serait pas possible/souhaitable de supprimer modd_lunit de PHYEX. On pourrait se contentner de recevoir le numero d'unité logique
+Regarder s'il ne serait pas possible/souhaitable de supprimer modd_lunit de PHYEX. On pourrait se contenter de recevoir le numero d'unité logique
Faire quelque chose de mesonh/micro/modd_blankn.f90: le déplacer dans common ou le supprimer
Nettoyage des répertoires aux nécessaire
+
+Initialiser dans AROME la variable ldiag_in_run de MODD_DIAG_IN_RUN pour pouvoir phaser le modd
diff --git a/src/arome/turb/mode_thermo_mono.F90 b/src/arome/aux/mode_thermo.F90
similarity index 98%
rename from src/arome/turb/mode_thermo_mono.F90
rename to src/arome/aux/mode_thermo.F90
index 8b85442eb8e537a9616ac40657d712c8c3d51aa7..86165d23b4df39acacf7c268d926aec1d7add298 100644
--- a/src/arome/turb/mode_thermo_mono.F90
+++ b/src/arome/aux/mode_thermo.F90
@@ -1,7 +1,10 @@
+!MNH_LIC Copyright 1994-2019 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.
+!-----------------------------------------------------------------
! #######################
MODULE MODE_THERMO
- USE PARKIND1, ONLY : JPRB
- USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! #######################
!
!!**** *MODE_THERMO_MONO* - module for routines SM_FOES,SM_PMR_HU
@@ -29,13 +32,23 @@
!! MODIFICATIONS
!! -------------
!! Original 28/08/94
+!! J.Escobar : 5/10/2018 : add FLUSH , for better logging in case of PB
+! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
!--------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
!-------------------------------------------------------------------------------
-!
+USE MODE_MSG
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: DQSAT, DQSATI, QSAT, QSATI, SM_FOES, SM_PMR_HU
+
INTERFACE SM_FOES
MODULE PROCEDURE SM_FOES_0D
MODULE PROCEDURE SM_FOES_1D
@@ -288,6 +301,8 @@ END FUNCTION SM_FOES_1D
!! Modification 16/03/95 remove the EPSILON function
!! Modification 15/09/97 (V. Masson) add solid and liquid water phases
!! in thetav computation
+!! Modification 22/01/2019 (P. Wautelet) use standard FLUSH statement
+!! instead of non standard intrinsics!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -347,7 +362,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_3D',0,ZHOOK_HANDLE)
ITERMAX = 10
IF (PRESENT(KITERMAX)) ITERMAX=KITERMAX
ZRDSRV = XRD /XRV
-ZEPS = 1.E-5
+ZEPS = XEPS_DT
!
ZRSLW(:,:,:)=0.
DO JRR=2,SIZE(PR,4)
@@ -379,8 +394,8 @@ IF ( ANY(ZDT > ZEPS) ) THEN
WRITE(ILUOUT,*) 'MR AT THIS MAXIMUM : ', PMR(IMAXLOC(1),IMAXLOC(2),IMAXLOC(3))
WRITE(ILUOUT,*) 'T AT THIS MAXIMUM : ', ZT(IMAXLOC(1),IMAXLOC(2),IMAXLOC(3))
WRITE(ILUOUT,*) 'JOB ABORTED '
- CALL ABORT
- STOP
+ FLUSH(unit=ILUOUT)
+ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'SM_PMR_HU_3D', 'failed to converge' )
END IF
!-------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_3D',1,ZHOOK_HANDLE)
@@ -522,8 +537,7 @@ IF (ANY(ZDT>ZEPS)) THEN
WRITE(ILUOUT,*) 'MR AT THIS MAXIMUM : ', PMR(IMAXLOC)
WRITE(ILUOUT,*) 'T AT THIS MAXIMUM : ', ZT(IMAXLOC)
WRITE(ILUOUT,*) 'JOB ABORTED '
- CALL ABORT
- STOP
+ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'SM_PMR_HU_1D', 'failed to converge' )
END IF
!-------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_1D',1,ZHOOK_HANDLE)
diff --git a/src/arome/turb/modd_dynn.F90 b/src/arome/dead_code/modd_dynn.F90
similarity index 100%
rename from src/arome/turb/modd_dynn.F90
rename to src/arome/dead_code/modd_dynn.F90
diff --git a/src/arome/turb/modn_turb.F90 b/src/arome/dead_code/modn_turb.F90
similarity index 100%
rename from src/arome/turb/modn_turb.F90
rename to src/arome/dead_code/modn_turb.F90
diff --git a/src/arome/gmkpack_ignored_files b/src/arome/gmkpack_ignored_files
index 0b47fe93506323efcd6859ed24a394ed7f57ff1e..193b04bd4624e14fe92f02563a87d775bf722f2d 100644
--- a/src/arome/gmkpack_ignored_files
+++ b/src/arome/gmkpack_ignored_files
@@ -164,3 +164,6 @@ phyex/micro/modd_nsv.F90
phyex/micro/modd_refaro.F90
phyex/micro/mode_fmbidon.F90
phyex/micro/mode_fmwritbidon.F90
+phyex/turb/modd_dynn.F90
+phyex/turb/mode_thermo_mono.F90
+phyex/turb/modn_turb.F90
diff --git a/src/arome/turb/ini_cturb.F90 b/src/arome/turb/ini_cturb.F90
index 06bbc545492660b517fea9f828a6abd13ee66937..1847d0fdaa16277c181230c08f4b43093856b0a7 100644
--- a/src/arome/turb/ini_cturb.F90
+++ b/src/arome/turb/ini_cturb.F90
@@ -1,7 +1,5 @@
! ######spl
SUBROUTINE INI_CTURB
- USE PARKIND1, ONLY : JPRB
- USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! ####################
!
!!**** *INI_CTURB* - routine to initialize the turbulence scheme
@@ -40,6 +38,8 @@
!! P.Jabouille 20/10/99 XCET=0.4
!! V.Masson 13/11/02 XALPSBL and XASBL
!! 05/06 Remove KEPS
+!! Q.Rodier 01/19
+!! Remove XASBL (not used)
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -48,8 +48,15 @@
USE MODD_CST
USE MODD_CTURB
!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
IMPLICIT NONE
!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('INI_CTURB',0,ZHOOK_HANDLE)
+!
! ---------------------------------------------------------------------------
!
! 1. SETTING THE NUMERICAL VALUES
@@ -57,8 +64,6 @@ IMPLICIT NONE
!
! 1.1 Constant for dissipation of Tke
!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('INI_CTURB',0,ZHOOK_HANDLE)
!
LHARAT=.FALSE.
!
@@ -67,6 +72,7 @@ XCED = 0.85
! Redelsperger-Sommeria (1981) = 0.70
! Schmidt-Schumann (1989) = 0.845
! Cheng-Canuto-Howard (2002) = 0.845
+! Rodier, Masson, Couvreux, Paci (2017) = 0.34
!
!
! 1.2 Constant for wind pressure-correlations
@@ -191,14 +197,6 @@ XCPR3= XCPR2 ! used only for the Schmidt number for scalar variables
XCPR4= XCPR2
XCPR5= XCPR2
!
-! 2.4 Value related to the TKE universal function within SBL
-!
-!
-XASBL = 0.5*( XALPSBL**(3./2.)*XKARMAN*XCED + XKARMAN/SQRT(XALPSBL)/XCMFS )
-! Redelsperger et al 2001
-!
-!
-!
! 3. MINIMUM VALUES
! --------------
!
diff --git a/src/arome/turb/modd_diag_in_run.F90 b/src/arome/turb/modd_diag_in_run.F90
index 48aa658f406545abf0fc72758a235ba5db5f9f16..ca9ce53a698a3a458f5ac777353e5d23518faa78 100644
--- a/src/arome/turb/modd_diag_in_run.F90
+++ b/src/arome/turb/modd_diag_in_run.F90
@@ -1,20 +1,30 @@
-! ######spl
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
MODULE MODD_DIAG_IN_RUN
+! Modifications
+!! 02/2018 Q.Libois ECRAD
+!! Bielli S. 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
!
!* stores instantaneous diagnostic arrays for the current time-step
!
IMPLICIT NONE
-SAVE
-
LOGICAL :: LDIAG_IN_RUN=.FALSE. ! flag for diagnostics
!
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_RN ! net radiation
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_H ! sensible heat flux
-REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LE ! latent heat flux
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LE ! Total latent heat flux
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LEI ! Solid latent heat flux
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_GFLUX ! ground flux
-REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LW ! incoming longwave at the surface
-REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SW ! incoming Shortwave at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LWD ! incoming longwave at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_LWU ! outcoming longwave at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SWD ! incoming Shortwave at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SWU ! outcoming Shortwave at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SWDIR ! incoming Shortwave direct at the surface
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SWDIFF! incoming Shortwave diffuse at the surface
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_T2M ! temperature at 2m
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_Q2M ! humidity at 2m
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_HU2M ! relative humidity at 2m
@@ -23,4 +33,6 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_MER10M! meridian wind at 10m
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_DSTAOD! dust aerosol optical depth
REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SFCO2 ! CO2 Surface flux
REAL, DIMENSION(:,:,:), ALLOCATABLE :: XCURRENT_TKE_DISS ! Tke dissipation rate
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_SLTAOD ! Salt aerosol optical depth
+REAL, DIMENSION(:,:), ALLOCATABLE :: XCURRENT_ZWS ! Significant height of waves
END MODULE MODD_DIAG_IN_RUN
diff --git a/src/arome/turb/modi_compute_function_thermo_mf.F90 b/src/arome/turb/modi_compute_function_thermo_mf.F90
deleted file mode 100644
index 7f68ce56fce49953f795cb474da461a774978b96..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_compute_function_thermo_mf.F90
+++ /dev/null
@@ -1,33 +0,0 @@
-! ######spl
- MODULE MODI_COMPUTE_FUNCTION_THERMO_MF
-! ######################################
-!
-INTERFACE
-
-! #################################################################
- SUBROUTINE COMPUTE_FUNCTION_THERMO_MF( KRR,KRRL,KRRI, &
- PTH, PR, PEXN, PFRAC_ICE, PPABS, &
- PT, PAMOIST,PATHETA )
-! #################################################################
-
-!* 1.1 Declaration of Arguments
-!
-
-INTEGER, INTENT(IN) :: KRR ! number of moist var.
-INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
-INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
-
-REAL, DIMENSION(:,:), INTENT(IN) :: PTH ! theta
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PR ! water species
-REAL, DIMENSION(:,:) , INTENT(IN) :: PPABS,PEXN ! pressure, Exner funct.
-REAL, DIMENSION(:,:) , INTENT(IN) :: PFRAC_ICE ! ice fraction
-
-REAL, DIMENSION(:,:), INTENT(OUT) :: PT ! temperature
-
-REAL, DIMENSION(:,:), INTENT(OUT) :: PAMOIST,PATHETA
-!
-END SUBROUTINE COMPUTE_FUNCTION_THERMO_MF
-
-END INTERFACE
-!
-END MODULE MODI_COMPUTE_FUNCTION_THERMO_MF
diff --git a/src/arome/turb/mode_sbl.F90 b/src/common/turb/mode_sbl.F90
similarity index 90%
rename from src/arome/turb/mode_sbl.F90
rename to src/common/turb/mode_sbl.F90
index bb24bb31e9c8bf8989f81a43ec06a3ac81291604..ef8e3ac6681eebf28d3df0b0351a6bd473149893 100644
--- a/src/arome/turb/mode_sbl.F90
+++ b/src/common/turb/mode_sbl.F90
@@ -1,7 +1,10 @@
-! ######spl
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+! ###############
MODULE MODE_SBL
- USE PARKIND1, ONLY : JPRB
- USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! ###############
!
!!**** *MODE_SBL * - contains Surface Boundary Layer characteristics functions
@@ -36,6 +39,8 @@
!! V. Masson 06/11/02 optimization and add Businger fonction for TKE
!! V. Masson 01/01/03 use PAULSON_PSIM function
!-----------------------------------------------------------------------------
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!
!* 0. DECLARATIONS
!
@@ -87,7 +92,7 @@ FUNCTION BUSINGER_PHIM_3D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIM_3D(:,:,:) = 1. + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_3D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_3D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIM_3D
!
!-------------------------------------------------------------------------------
@@ -103,7 +108,7 @@ FUNCTION BUSINGER_PHIM_2D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIM_2D(:,:) = 1. + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_2D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_2D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIM_2D
!
!-------------------------------------------------------------------------------
@@ -119,7 +124,7 @@ FUNCTION BUSINGER_PHIM_1D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIM_1D(:) = 1. + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_1D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_1D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIM_1D
!
!-------------------------------------------------------------------------------
@@ -135,7 +140,7 @@ FUNCTION BUSINGER_PHIM_0D(PZ_O_LMO)
ELSE
BUSINGER_PHIM_0D = 1. + 4.7 * PZ_O_LMO
END IF
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_0D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM_0D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIM_0D
!
!-------------------------------------------------------------------------------
@@ -153,7 +158,7 @@ FUNCTION BUSINGER_PHIH_3D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIH_3D(:,:,:) = 0.74 + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_3D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_3D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIH_3D
!
!-------------------------------------------------------------------------------
@@ -169,7 +174,7 @@ FUNCTION BUSINGER_PHIH_2D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIH_2D(:,:) = 0.74 + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_2D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_2D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIH_2D
!
!-------------------------------------------------------------------------------
@@ -185,7 +190,7 @@ FUNCTION BUSINGER_PHIH_1D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIH_1D(:) = 0.74 + 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_1D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_1D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIH_1D
!
!-------------------------------------------------------------------------------
@@ -201,7 +206,7 @@ FUNCTION BUSINGER_PHIH_0D(PZ_O_LMO)
ELSE
BUSINGER_PHIH_0D = 0.74 + 4.7 * PZ_O_LMO
END IF
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_0D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH_0D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIH_0D
!
!-------------------------------------------------------------------------------
@@ -221,7 +226,7 @@ FUNCTION BUSINGER_PHIE_3D(PZ_O_LMO)
ELSEWHERE
BUSINGER_PHIE_3D(:,:,:) = 1./(1. + 4.7 * PZ_O_LMO)**2
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIE_3D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIE_3D',1,ZHOOK_HANDLE)
END FUNCTION BUSINGER_PHIE_3D
!
!-------------------------------------------------------------------------------
@@ -243,7 +248,7 @@ FUNCTION PAULSON_PSIM_2D(PZ_O_LMO)
ELSEWHERE
PAULSON_PSIM_2D(:,:) = - 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_2D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_2D',1,ZHOOK_HANDLE)
END FUNCTION PAULSON_PSIM_2D
!
!-------------------------------------------------------------------------------
@@ -264,7 +269,7 @@ FUNCTION PAULSON_PSIM_1D(PZ_O_LMO)
ELSEWHERE
PAULSON_PSIM_1D(:) = - 4.7 * PZ_O_LMO
END WHERE
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_1D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_1D',1,ZHOOK_HANDLE)
END FUNCTION PAULSON_PSIM_1D
!
!-------------------------------------------------------------------------------
@@ -285,7 +290,7 @@ FUNCTION PAULSON_PSIM_0D(PZ_O_LMO)
ELSE
PAULSON_PSIM_0D = - 4.7 * PZ_O_LMO
END IF
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_0D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:PAULSON_PSIM_0D',1,ZHOOK_HANDLE)
END FUNCTION PAULSON_PSIM_0D
!
!-------------------------------------------------------------------------------
@@ -317,7 +322,7 @@ FUNCTION LMO_2D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
LMO_2D(:,:) = - MAX(PUSTAR(:,:),1.E-6)**3 &
/ ( XKARMAN * XG / ZTHETAV(:,:) *ZQ0(:,:) )
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_2D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_2D',1,ZHOOK_HANDLE)
END FUNCTION LMO_2D
!
!-------------------------------------------------------------------------------
@@ -347,7 +352,7 @@ FUNCTION LMO_1D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
LMO_1D(:) = - MAX(PUSTAR(:),1.E-6)**3 &
/ ( XKARMAN * ( XG / ZTHETAV(:) * PSFTH(:) &
+ XG * ZEPS * PSFRV(:) ) )
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_1D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_1D',1,ZHOOK_HANDLE)
END FUNCTION LMO_1D
!
!-------------------------------------------------------------------------------
@@ -378,7 +383,7 @@ FUNCTION LMO_0D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
LMO_0D = - MAX(PUSTAR,1.E-6)**3 &
/ ( XKARMAN * ( XG / ZTHETAV * PSFTH &
+ XG * ZEPS * PSFRV ) )
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_0D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO_0D',1,ZHOOK_HANDLE)
END FUNCTION LMO_0D
!
!-------------------------------------------------------------------------------
@@ -421,7 +426,7 @@ FUNCTION USTAR_2D(PU,PV,PZ,PZ0,PLMO)
* XKARMAN / LOG(PZ(:,:)/PZ0(:,:))
END WHERE
!
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_2D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_2D',1,ZHOOK_HANDLE)
END FUNCTION USTAR_2D
!
!-------------------------------------------------------------------------------
@@ -463,7 +468,7 @@ FUNCTION USTAR_1D(PU,PV,PZ,PZ0,PLMO)
* XKARMAN / LOG(PZ(:)/PZ0(:))
END WHERE
!
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_1D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_1D',1,ZHOOK_HANDLE)
END FUNCTION USTAR_1D
!
!-------------------------------------------------------------------------------
@@ -495,7 +500,7 @@ FUNCTION USTAR_0D(PU,PV,PZ,PZ0,PLMO)
USTAR_0D = SQRT( PU**2+PV**2 ) &
* XKARMAN / LOG(PZ/PZ0)
-IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_0D',1,ZHOOK_HANDLE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:USTAR_0D',1,ZHOOK_HANDLE)
END FUNCTION USTAR_0D
!
!-------------------------------------------------------------------------------
diff --git a/src/arome/turb/modi_les_mean_subgrid.F90 b/src/common/turb/modi_les_mean_subgrid.F90
similarity index 100%
rename from src/arome/turb/modi_les_mean_subgrid.F90
rename to src/common/turb/modi_les_mean_subgrid.F90
diff --git a/src/mesonh/aux/mode_thermo.f90 b/src/mesonh/aux/mode_thermo.f90
new file mode 100644
index 0000000000000000000000000000000000000000..935ffebd0606833ea3979137886485da05f041bd
--- /dev/null
+++ b/src/mesonh/aux/mode_thermo.f90
@@ -0,0 +1,2515 @@
+!MNH_LIC Copyright 1994-2019 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.
+!-----------------------------------------------------------------
+! #######################
+ MODULE MODE_THERMO
+! #######################
+!
+!!**** *MODE_THERMO_MONO* - module for routines SM_FOES,SM_PMR_HU
+!!
+!! PURPOSE
+!! -------
+! The purpose of this executive module is to package
+! the routine SM_FOES, SM_PMR_HU without use of comlib parallel routine
+!
+!
+!
+!!
+!!** IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! J.Escobar : 5/10/2018 : add FLUSH , for better logging in case of PB
+! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!--------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!-------------------------------------------------------------------------------
+USE MODE_MSG
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: DQSAT, DQSATI, QSAT, QSATI, SM_FOES, SM_PMR_HU
+
+INTERFACE SM_FOES
+ MODULE PROCEDURE SM_FOES_0D
+ MODULE PROCEDURE SM_FOES_1D
+ MODULE PROCEDURE SM_FOES_2D
+ MODULE PROCEDURE SM_FOES_2D_MASK
+ MODULE PROCEDURE SM_FOES_3D
+END INTERFACE
+INTERFACE QSAT
+ MODULE PROCEDURE QSATW_3D
+ MODULE PROCEDURE QSATW_2D
+ MODULE PROCEDURE QSATW_2D_MASK
+ MODULE PROCEDURE QSATW_1D
+ MODULE PROCEDURE QSATW_0D
+END INTERFACE
+INTERFACE DQSAT
+ MODULE PROCEDURE DQSATW_O_DT_2D_MASK
+ MODULE PROCEDURE DQSATW_O_DT_1D
+ MODULE PROCEDURE DQSATW_O_DT_3D
+END INTERFACE
+INTERFACE QSATI
+ MODULE PROCEDURE QSATI_3D
+ MODULE PROCEDURE QSATI_2D
+ MODULE PROCEDURE QSATI_2D_MASK
+ MODULE PROCEDURE QSATI_1D
+ MODULE PROCEDURE QSATI_0D
+END INTERFACE
+INTERFACE DQSATI
+ MODULE PROCEDURE DQSATI_O_DT_2D_MASK
+ MODULE PROCEDURE DQSATI_O_DT_1D
+ MODULE PROCEDURE DQSATI_O_DT_3D
+END INTERFACE
+INTERFACE SM_PMR_HU
+ MODULE PROCEDURE SM_PMR_HU_1D
+ MODULE PROCEDURE SM_PMR_HU_3D
+END INTERFACE
+CONTAINS
+!-------------------------------------------------------------------------------
+! ####################################
+ FUNCTION SM_FOES_3D(PT) RESULT(PFOES)
+! ####################################
+!
+!!**** *SM_FOES_3D * - function to compute saturation vapor pressure from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: PFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!* 0.2 Declarations of local variables
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_3D',0,ZHOOK_HANDLE)
+PFOES(:,:,:) = EXP( XALPW - XBETAW/PT(:,:,:) - XGAMW*LOG(PT(:,:,:)) )
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_3D',1,ZHOOK_HANDLE)
+END FUNCTION SM_FOES_3D
+! ####################################
+ FUNCTION SM_FOES_1D(PT) RESULT(PFOES)
+! ####################################
+!
+!!**** *SM_FOES_1D * - function to compute saturation vapor pressure from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(SIZE(PT)) :: PFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!* 0.2 Declarations of local variables
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_1D',0,ZHOOK_HANDLE)
+PFOES(:) = EXP( XALPW - XBETAW/PT(:) - XGAMW*LOG(PT(:)) )
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_1D',1,ZHOOK_HANDLE)
+END FUNCTION SM_FOES_1D
+!-------------------------------------------------------------------------------
+! ####################################################
+ FUNCTION SM_PMR_HU_3D(PP,PTV,PHU,PR,KITERMAX) RESULT(PMR)
+! ####################################################
+!
+!!**** *SM_PMR_HU_3D * - function to compute vapor mixing ratio
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the vapor mixing ratio
+! from pressure, virtual temperature and relative humidity
+!
+!
+!!** METHOD
+!! ------
+!! Given Pressure (PP), Virtual temperature (PTV) and Relative
+!! humidity (PHU), the vapor mixing ratio is computed by iterating
+!! the following procedure :
+!! T ----> es(T)
+!! es(T) ,HU ----> es(Td)
+!! es(Td), P ----> r
+!! r , Tv ----> T
+!!
+!! at the beginning T=Tv
+!!
+!! EXTERNAL
+!! --------
+!! SM_FOES : to compute saturation vapor pressure
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XRV : gas constant for vapor
+!! XRD : gas constant for dry air
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 30/08/94
+!! Modification 16/03/95 remove the EPSILON function
+!! Modification 15/09/97 (V. Masson) add solid and liquid water phases
+!! in thetav computation
+!! Modification 22/01/2019 (P. Wautelet) use standard FLUSH statement
+!! instead of non standard intrinsics!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_LUNIT_n, ONLY: TLUOUT
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTV ! Virtual Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PHU ! Relative humidity
+ ! (percent)
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR ! vapor, liquid and
+ ! solid water mixing
+ ! ratio
+
+INTEGER, INTENT(IN), OPTIONAL :: KITERMAX ! maximum number
+ ! of iterations
+ ! (default 10)
+!
+REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2),SIZE(PP,3)) :: PMR ! vapor mixing ratio
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2),SIZE(PP,3)) :: ZT ! temperature
+REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2),SIZE(PP,3)) :: ZDT ! increment of
+ ! temperature between two iterations
+REAL :: ZRDSRV ! Rd/Rv
+REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2),SIZE(PP,3)) :: ZESTD ! es(Td)
+REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2),SIZE(PP,3)) :: ZRSLW ! total solid and liquid water mixing ratio
+INTEGER :: ITERMAX ! Maximum number
+ ! of iteration
+INTEGER :: ITER ! iteration number of
+REAL :: ZEPS ! a small number
+INTEGER, DIMENSION(3) :: IMAXLOC ! localisation of
+ ! a maximum
+INTEGER :: ILUOUT
+ ! logical unit for
+ ! output-listing
+ ! and error code
+INTEGER :: JRR ! loop counter
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE VAPOR MIXING RATIO
+! --------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_3D',0,ZHOOK_HANDLE)
+ITERMAX = 10
+IF (PRESENT(KITERMAX)) ITERMAX=KITERMAX
+ZRDSRV = XRD /XRV
+ZEPS = XEPS_DT
+!
+ZRSLW(:,:,:)=0.
+DO JRR=2,SIZE(PR,4)
+ ZRSLW(:,:,:)=ZRSLW(:,:,:)+PR(:,:,:,JRR)
+END DO
+!
+ZT(:,:,:) = PTV(:,:,:)
+DO ITER=1,ITERMAX
+ ZESTD(:,:,:) = PHU(:,:,:) * SM_FOES(ZT(:,:,:)) * 0.01
+ PMR (:,:,:) = ZRDSRV * ZESTD(:,:,:) /(PP(:,:,:) - ZESTD(:,:,:))
+ ZDT(:,:,:) = ZT(:,:,:)
+ ZT(:,:,:) = PTV(:,:,:) * (1.+PMR(:,:,:)+ZRSLW(:,:,:)) / (1.+ PMR(:,:,:)/ZRDSRV)
+ ZDT(:,:,:) = ABS(ZDT(:,:,:) - ZT(:,:,:))
+END DO
+!-------------------------------------------------------------------------------
+!
+!* 2. NO CONVERGENCE
+! --------------
+!
+IF ( ANY(ZDT > ZEPS) ) THEN
+ ILUOUT = TLUOUT%NLU
+ WRITE(ILUOUT,*) 'ERROR IN FUNCTION SM_PMR_HU (module MODE_THERMO)'
+ WRITE(ILUOUT,*) 'FUNCTION FAILS TO CONVERGE AFTER ', ITERMAX,' ITERATIONS'
+ WRITE(ILUOUT,*) 'EPS = ' , ZEPS
+ IMAXLOC(:) = MAXLOC(ZDT)
+ WRITE(ILUOUT,*) 'MAXIMUM RESIDUAL DT :', MAXVAL(ZDT)
+! WRITE(ILUOUT,*) 'LOCATION OF THIS MAXIMUM I=',IMAXLOC(1),' J=',IMAXLOC(2), &
+! ' K=',IMAXLOC(3)
+ WRITE(ILUOUT,*) 'MR AT THIS MAXIMUM : ', PMR(IMAXLOC(1),IMAXLOC(2),IMAXLOC(3))
+ WRITE(ILUOUT,*) 'T AT THIS MAXIMUM : ', ZT(IMAXLOC(1),IMAXLOC(2),IMAXLOC(3))
+ WRITE(ILUOUT,*) 'JOB ABORTED '
+ FLUSH(unit=ILUOUT)
+ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'SM_PMR_HU_3D', 'failed to converge' )
+END IF
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_3D',1,ZHOOK_HANDLE)
+END FUNCTION SM_PMR_HU_3D
+! ################################################################
+ FUNCTION SM_PMR_HU_1D(PP,PTV,PHU,PR,KITERMAX) RESULT(PMR)
+! ################################################################
+!
+!!**** *SM_PMR_HU_1D * - function to compute vapor mixing ratio
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the vapor mixing ratio
+! from pressure, virtual temperature and relative humidity
+!
+!
+!!** METHOD
+!! ------
+!! Given Pressure (PP), Virtual temperature (PTV) and Relative
+!! humidity (PHU), the vapor mixing ratio is computed by iterating
+!! the following procedure :
+!! T ----> es(T)
+!! es(T) ,HU ----> es(Td)
+!! es(Td), P ----> r
+!! r , Tv ----> T
+!!
+!! at the beginning T=Tv
+!!
+!! EXTERNAL
+!! --------
+!! SM_FOES : to compute saturation vapor pressure
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XRV : gas constant for vapor
+!! XRD : gas constant for dry air
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 30/08/94
+!! Modification 16/03/95 remove the EPSILON function
+!! Modification 15/09/97 (V. Masson) add solid and liquid water phases
+!! in thetav computation
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_LUNIT_n, ONLY: TLUOUT
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:), INTENT(IN) :: PTV ! Virtual Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:), INTENT(IN) :: PHU ! Relative humidity
+ ! (percent)
+REAL, DIMENSION(:,:), INTENT(IN) :: PR ! vapor, liquid and solid
+ ! water mixing ratio
+INTEGER, INTENT(IN), OPTIONAL :: KITERMAX ! maximum number
+ ! of iterations
+ ! (default 10)
+!
+REAL, DIMENSION(SIZE(PP)) :: PMR ! vapor mixing ratio
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PP)) :: ZT ! temperature
+REAL, DIMENSION(SIZE(PP)) :: ZDT ! increment of
+ ! temperature between two iterations
+REAL :: ZRDSRV ! Rd/Rv
+REAL, DIMENSION(SIZE(PP)) :: ZESTD ! es(Td)
+REAL, DIMENSION(SIZE(PP)) :: ZRSLW ! total solid and liquid water mixing ratio
+INTEGER :: ITERMAX ! Maximum number
+ ! of iteration
+INTEGER :: ITER ! iteration number of
+REAL :: ZEPS ! a small number
+INTEGER,DIMENSION(1) :: IMAXLOC ! localisation of
+ ! a maximum
+INTEGER :: ILUOUT,IRESP
+ ! logical unit for
+ ! output-listing
+ ! and error code
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE VAPOR MIXING RATIO
+! --------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_1D',0,ZHOOK_HANDLE)
+ITERMAX = 10
+IF (PRESENT(KITERMAX)) ITERMAX=KITERMAX
+ZRDSRV = XRD /XRV
+ZEPS = 1.E-5
+!
+IF (SIZE(PR,2)>1) THEN
+ ZRSLW(:)=SUM(PR(:,2:),DIM=2)
+ELSE
+ ZRSLW(:)=0.
+END IF
+!
+ZT(:) = PTV(:)
+DO ITER=1,ITERMAX
+ ZESTD(:) = PHU(:) * SM_FOES(ZT(:)) * 0.01
+ PMR (:) = ZRDSRV * ZESTD(:) /(PP(:) - ZESTD(:))
+ ZDT(:) = ZT(:)
+ ZT(:) = PTV(:) * (1.+PMR(:)+ZRSLW(:)) / (1.+ PMR(:)/ZRDSRV)
+ ZDT(:) = ABS(ZDT(:) - ZT(:))
+END DO
+!-------------------------------------------------------------------------------
+!
+!* 2. NO CONVERGENCE
+! --------------
+!
+IF (ANY(ZDT>ZEPS)) THEN
+ ILUOUT = TLUOUT%NLU
+ WRITE(ILUOUT,*) 'ERROR IN FUNCTION SM_PMR_HU (module MODE_THERMO)'
+ WRITE(ILUOUT,*) 'FUNCTION FAILS TO CONVERGE AFTER ', ITERMAX,' ITERATIONS'
+ WRITE(ILUOUT,*) 'EPS = ' , ZEPS
+ IMAXLOC = MAXLOC(ZDT)
+ WRITE(ILUOUT,*) 'MAXIMUM RESIDUAL DT :', MAXVAL(ZDT)
+ WRITE(ILUOUT,*) 'MR AT THIS MAXIMUM : ', PMR(IMAXLOC)
+ WRITE(ILUOUT,*) 'T AT THIS MAXIMUM : ', ZT(IMAXLOC)
+ WRITE(ILUOUT,*) 'JOB ABORTED '
+ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'SM_PMR_HU_1D', 'failed to converge' )
+END IF
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_PMR_HU_1D',1,ZHOOK_HANDLE)
+END FUNCTION SM_PMR_HU_1D
+! ####################################
+ FUNCTION SM_FOES_0D(PT) RESULT(PFOES)
+! ####################################
+!
+!!**** *SM_FOES_0D * - function to compute saturation vapor pressure from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!! 24/12/97 (V. Masson) version for 0D arrays
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL :: PFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!* 0.2 Declarations of local variables
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_0D',0,ZHOOK_HANDLE)
+PFOES = EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT) )
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_0D',1,ZHOOK_HANDLE)
+END FUNCTION SM_FOES_0D
+!
+!-------------------------------------------------------------------------------
+! ####################################
+ FUNCTION SM_FOES_2D(PT) RESULT(PFOES)
+! ####################################
+!
+!!**** *SM_FOES_2D * - function to compute saturation vapor pressure from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!! 24/12/97 (V. Masson) version for 2D arrays
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!* 0.2 Declarations of local variables
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_2D',0,ZHOOK_HANDLE)
+PFOES(:,:) = EXP( XALPW - XBETAW/PT(:,:) - XGAMW*LOG(PT(:,:)) )
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_2D',1,ZHOOK_HANDLE)
+END FUNCTION SM_FOES_2D
+!
+!-------------------------------------------------------------------------------
+!
+! ################################################
+ FUNCTION SM_FOES_2D_MASK(OMASK,PT) RESULT(PFOES)
+! ################################################
+!
+!!**** *SM_FOES_2D * - function to compute saturation vapor pressure from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Ducrocq * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 28/08/94
+!! 24/12/97 (V. Masson) version for 2D arrays
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK ! Localization mask
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!* 0.2 Declarations of local variables
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_2D_MASK',0,ZHOOK_HANDLE)
+WHERE (OMASK(:,:))
+ PFOES(:,:) = EXP( XALPW - XBETAW/PT(:,:) - XGAMW*LOG(PT(:,:)) )
+END WHERE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:SM_FOES_2D_MASK',1,ZHOOK_HANDLE)
+END FUNCTION SM_FOES_2D_MASK
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATW_3D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_3D',0,ZHOOK_HANDLE)
+ZFOES(:,:,:) = MIN(EXP( XALPW - XBETAW/PT(:,:,:) - XGAMW*LOG(PT(:,:,:)) ), 0.99*PP(:,:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:,:,:) = XRD/XRV*ZFOES(:,:,:)/PP(:,:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:,:)/PP(:,:,:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_3D',1,ZHOOK_HANDLE)
+END FUNCTION QSATW_3D
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATW_2D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_2D',0,ZHOOK_HANDLE)
+ZFOES(:,:) = MIN(EXP( XALPW - XBETAW/PT(:,:) - XGAMW*LOG(PT(:,:)) ), 0.99*PP(:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_2D',1,ZHOOK_HANDLE)
+END FUNCTION QSATW_2D
+!
+!-------------------------------------------------------------------------------
+!
+! #################################################
+ FUNCTION QSATW_2D_MASK(OMASK,PT,PP) RESULT(PQSAT)
+! #################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK ! Localization mask
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_2D_MASK',0,ZHOOK_HANDLE)
+WHERE (OMASK(:,:))
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ ZFOES(:,:) = MIN(EXP( XALPW - XBETAW/PT(:,:) - XGAMW*LOG(PT(:,:)) ), 0.99*PP(:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+ PQSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
+ELSEWHERE
+!
+!* 3. BOGUS VALUE
+! -----------
+!
+ PQSAT(:,:) = 0.
+END WHERE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_2D_MASK',1,ZHOOK_HANDLE)
+END FUNCTION QSATW_2D_MASK
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATW_1D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_1D',0,ZHOOK_HANDLE)
+ZFOES(:) = MIN(EXP( XALPW - XBETAW/PT(:) - XGAMW*LOG(PT(:)) ), 0.99*PP(:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:) = XRD/XRV*ZFOES(:)/PP(:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_1D',1,ZHOOK_HANDLE)
+END FUNCTION QSATW_1D
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATW_0D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_0D',0,ZHOOK_HANDLE)
+ZFOES = MIN(EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT) ), 0.99*PP)
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT = XRD/XRV*ZFOES/PP / (1.+(XRD/XRV-1.)*ZFOES/PP)
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATW_0D',1,ZHOOK_HANDLE)
+END FUNCTION QSATW_0D
+!
+!-------------------------------------------------------------------------------
+!
+! ##############################################################
+ FUNCTION DQSATW_O_DT_2D_MASK(OMASK,PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK ! Localization mask
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:,:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATW_O_DT_2D_MASK',0,ZHOOK_HANDLE)
+WHERE (OMASK(:,:))
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ ZFOES(:,:) = PP(:,:) / (1.+XRD/XRV*(1./PQSAT(:,:)-1.))
+!
+!* 2. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+ PDQSAT(:,:) = PQSAT(:,:) / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:) ) &
+ * (XBETAW/PT(:,:)**2 - XGAMW/PT(:,:))
+ELSEWHERE
+!
+!* 3. BOGUS VALUE
+! -----------
+!
+ PDQSAT(:,:) = 0.
+END WHERE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATW_O_DT_2D_MASK',1,ZHOOK_HANDLE)
+END FUNCTION DQSATW_O_DT_2D_MASK
+!
+!-------------------------------------------------------------------------------
+! ##############################################################
+ FUNCTION DQSATW_O_DT_1D(PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATW_O_DT_1D',0,ZHOOK_HANDLE)
+ZFOES(:) = PP(:) / (1.+XRD/XRV*(1./PQSAT(:)-1.))
+!
+!* 2. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+PDQSAT(:) = PQSAT(:) / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:) ) &
+ * (XBETAW/PT(:)**2 - XGAMW/PT(:))
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATW_O_DT_1D',1,ZHOOK_HANDLE)
+END FUNCTION DQSATW_O_DT_1D
+!
+!-------------------------------------------------------------------------------
+! ##############################################################
+ FUNCTION DQSATW_O_DT_3D(PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ZFOES(:,:,:) = PP(:,:,:) / (1.+XRD/XRV*(1./PQSAT(:,:,:)-1.))
+!
+!* 2. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+PDQSAT(:,:,:) = PQSAT(:,:,:) / (1.+(XRD/XRV-1.)*ZFOES(:,:,:)/PP(:,:,:) ) &
+ * (XBETAW/PT(:,:,:)**2 - XGAMW/PT(:,:,:))
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DQSATW_O_DT_3D
+!
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+!
+! ##############################################################
+ FUNCTION DQSATI_O_DT_2D_MASK(OMASK,PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature (with respect to ice)
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK ! Localization mask
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:,:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATI_O_DT_2D_MASK',0,ZHOOK_HANDLE)
+WHERE (OMASK(:,:))
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ ZFOES(:,:) = PP(:,:) / (1.+XRD/XRV*(1./PQSAT(:,:)-1.))
+!
+!* 3. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+ PDQSAT(:,:) = PQSAT(:,:) / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:) ) &
+ * (XBETAI/PT(:,:)**2 - XGAMI/PT(:,:))
+ELSEWHERE
+!
+!* 3. BOGUS VALUE
+! -----------
+!
+ PDQSAT(:,:) = 0.
+END WHERE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATI_O_DT_2D_MASK',1,ZHOOK_HANDLE)
+END FUNCTION DQSATI_O_DT_2D_MASK
+!
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! ##############################################################
+ FUNCTION DQSATI_O_DT_1D(PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature (with respect to ice)
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATI_O_DT_1D',0,ZHOOK_HANDLE)
+ZFOES(:) = PP(:) / (1.+XRD/XRV*(1./PQSAT(:)-1.))
+!
+!* 3. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+PDQSAT(:) = PQSAT(:) / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:) ) &
+ * (XBETAI/PT(:)**2 - XGAMI/PT(:))
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:DQSATI_O_DT_1D',1,ZHOOK_HANDLE)
+END FUNCTION DQSATI_O_DT_1D
+!
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! ##############################################################
+ FUNCTION DQSATI_O_DT_3D(PT,PP,PQSAT) RESULT(PDQSAT)
+! ##############################################################
+!
+!!**** *QSATW * - function to compute saturation vapor humidity from
+!! temperature (with respect to ice)
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!! Finally, dqsat / dT (T) is computed.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPW : Constant for saturation vapor pressure function
+!! XBETAW : Constant for saturation vapor pressure function
+!! XGAMW : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: PDQSAT ! derivative according
+ ! to temperature of
+ ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg))
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ZFOES(:,:,:) = PP(:,:,:) / (1.+XRD/XRV*(1./PQSAT(:,:,:)-1.))
+!
+!* 3. DERIVATION ACCORDING TO TEMPERATURE
+! -----------------------------------
+!
+PDQSAT(:,:,:) = PQSAT(:,:,:) / (1.+(XRD/XRV-1.)*ZFOES(:,:,:)/PP(:,:,:) ) &
+ * (XBETAI/PT(:,:,:)**2 - XGAMI/PT(:,:,:))
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION DQSATI_O_DT_3D
+!
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATI_3D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATI * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPI) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAI) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMI) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPI : Constant for saturation vapor pressure function
+!! XBETAI : Constant for saturation vapor pressure function
+!! XGAMI : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2),SIZE(PT,3)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_3D',0,ZHOOK_HANDLE)
+ZFOES(:,:,:) = MIN(EXP( XALPI - XBETAI/PT(:,:,:) - XGAMI*LOG(PT(:,:,:)) ), 0.99*PP(:,:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:,:,:) = XRD/XRV*ZFOES(:,:,:)/PP(:,:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:,:)/PP(:,:,:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_3D',1,ZHOOK_HANDLE)
+END FUNCTION QSATI_3D
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATI_2D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATI * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPI) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAI) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMI) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPI : Constant for saturation vapor pressure function
+!! XBETAI : Constant for saturation vapor pressure function
+!! XGAMI : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_2D',0,ZHOOK_HANDLE)
+ZFOES(:,:) = MIN(EXP( XALPI - XBETAI/PT(:,:) - XGAMI*LOG(PT(:,:)) ), 0.99*PP(:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_2D',1,ZHOOK_HANDLE)
+END FUNCTION QSATI_2D
+!
+!-------------------------------------------------------------------------------
+!
+! #################################################
+ FUNCTION QSATI_2D_MASK(OMASK,PT,PP) RESULT(PQSAT)
+! #################################################
+!
+!!**** *QSATI * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPI) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAI) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMI) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPI : Constant for saturation vapor pressure function
+!! XBETAI : Constant for saturation vapor pressure function
+!! XGAMI : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK ! Localization mask
+REAL, DIMENSION(:,:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1),SIZE(PT,2)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_2D_MASK',0,ZHOOK_HANDLE)
+WHERE (OMASK(:,:))
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+ ZFOES(:,:) = MIN(EXP( XALPI - XBETAI/PT(:,:) - XGAMI*LOG(PT(:,:)) ), 0.99*PP(:,:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+ PQSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
+ELSEWHERE
+!
+!* 3. BOGUS VALUE
+! -----------
+!
+ PQSAT(:,:) = 0.
+END WHERE
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_2D_MASK',1,ZHOOK_HANDLE)
+END FUNCTION QSATI_2D_MASK
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATI_1D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATI * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPI) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAI) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMI) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPI : Constant for saturation vapor pressure function
+!! XBETAI : Constant for saturation vapor pressure function
+!! XGAMI : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, DIMENSION(:), INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, DIMENSION(:), INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL, DIMENSION(SIZE(PT,1)) :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT,1)) :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_1D',0,ZHOOK_HANDLE)
+ZFOES(:) = MIN(EXP( XALPI - XBETAI/PT(:) - XGAMI*LOG(PT(:)) ), 0.99*PP(:))
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT(:) = XRD/XRV*ZFOES(:)/PP(:) &
+ / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:))
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_1D',1,ZHOOK_HANDLE)
+END FUNCTION QSATI_1D
+!
+!-------------------------------------------------------------------------------
+!
+! ######################################
+ FUNCTION QSATI_0D(PT,PP) RESULT(PQSAT)
+! ######################################
+!
+!!**** *QSATI * - function to compute saturation vapor humidity from
+!! temperature
+!!
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the saturation vapor
+! pressure from temperature
+!
+!
+!!** METHOD
+!! ------
+!! Given temperature T (PT), the saturation vapor pressure es(T)
+!! (FOES(PT)) is computed by integration of the Clapeyron equation
+!! from the triple point temperature Tt (XTT) and the saturation vapor
+!! pressure of the triple point es(Tt) (XESTT), i.e
+!!
+!! es(T)= EXP( alphaw - betaw /T - gammaw Log(T) )
+!!
+!! with :
+!! alphaw (XALPI) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt)
+!! betaw (XBETAI) = Lv(Tt)/Rv + gammaw Tt
+!! gammaw (XGAMI) = (Cl -Cpv) /Rv
+!!
+!! Then, the specific humidity at saturation is deduced.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : comtains physical constants
+!! XALPI : Constant for saturation vapor pressure function
+!! XBETAI : Constant for saturation vapor pressure function
+!! XGAMI : Constant for saturation vapor pressure function
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 21/09/98
+!! S. Riette april 2011 : protection in high statosphere where ZFOES > PP
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+!
+!
+REAL, INTENT(IN) :: PT ! Temperature
+ ! (Kelvin)
+REAL, INTENT(IN) :: PP ! Pressure
+ ! (Pa)
+REAL :: PQSAT ! saturation vapor
+ ! specific humidity
+ ! with respect to
+ ! water (kg/kg)
+!
+!* 0.2 Declarations of local variables
+!
+REAL :: ZFOES ! saturation vapor
+ ! pressure
+ ! (Pascal)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE SATURATION VAPOR PRESSURE
+! ---------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_0D',0,ZHOOK_HANDLE)
+ZFOES = MIN(EXP( XALPI - XBETAI/PT - XGAMI*LOG(PT) ), 0.99*PP)
+!
+!* 2. COMPUTE SATURATION HUMIDITY
+! ---------------------------
+!
+PQSAT = XRD/XRV*ZFOES/PP / (1.+(XRD/XRV-1.)*ZFOES/PP)
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('MODE_THERMO:QSATI_0D',1,ZHOOK_HANDLE)
+END FUNCTION QSATI_0D
+!
+!-------------------------------------------------------------------------------
+END MODULE MODE_THERMO
diff --git a/src/mesonh/turb/ini_cturb.f90 b/src/mesonh/turb/ini_cturb.f90
index 245dfa0632b533d2855e56b9fc3459a9e51a48cd..4da0d6452f856d60d7f795669619c3eca4dcc6ba 100644
--- a/src/mesonh/turb/ini_cturb.f90
+++ b/src/mesonh/turb/ini_cturb.f90
@@ -73,8 +73,15 @@ END MODULE MODI_INI_CTURB
USE MODD_CST
USE MODD_CTURB
!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
IMPLICIT NONE
!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('INI_CTURB',0,ZHOOK_HANDLE)
+!
! ---------------------------------------------------------------------------
!
! 1. SETTING THE NUMERICAL VALUES
@@ -85,7 +92,6 @@ IMPLICIT NONE
!XCED is now replaced by XCEDIS
!XCED = 0.70
!XCED = 0.84
-!
! Redelsperger-Sommeria (1981) = 0.70
! Schmidt-Schumann (1989) = 0.845
! Cheng-Canuto-Howard (2002) = 0.845
@@ -251,4 +257,5 @@ XSBL_O_BL = 0.05 ! SBL height / BL height ratio
XFTOP_O_FSURF = 0.05 ! Fraction of surface (heat or momentum) flux used to define top of BL
!
!
+IF (LHOOK) CALL DR_HOOK('INI_CTURB',1,ZHOOK_HANDLE)
END SUBROUTINE INI_CTURB
diff --git a/src/mesonh/turb/modd_diag_in_run.f90 b/src/mesonh/turb/modd_diag_in_run.f90
index b7bba80d0c045a7787cf64d952de4b44c6a2961f..6f9829570ec8a9f75685491317b5adc86b18e623 100644
--- a/src/mesonh/turb/modd_diag_in_run.f90
+++ b/src/mesonh/turb/modd_diag_in_run.f90
@@ -3,11 +3,6 @@
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 modd 2006/10/24 10:07:40
-!-----------------------------------------------------------------
MODULE MODD_DIAG_IN_RUN
! Modifications
!! 02/2018 Q.Libois ECRAD
diff --git a/src/mesonh/turb/mode_sbl.f90 b/src/mesonh/turb/mode_sbl.f90
deleted file mode 100644
index 1c5e1da7f4600a7d6fe0a5a228b4708bab55c526..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/mode_sbl.f90
+++ /dev/null
@@ -1,457 +0,0 @@
-!MNH_LIC Copyright 1994-2014 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.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 mode 2006/05/18 13:07:25
-!-----------------------------------------------------------------
-! ###############
- MODULE MODE_SBL
-! ###############
-!
-!!**** *MODE_SBL * - contains Surface Boundary Layer characteristics functions
-!!
-!! PURPOSE
-!! -------
-!
-!!** METHOD
-!! ------
-!!
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!!
-!! Businger et al 1971, Wyngaard and Cote 1974
-!!
-!!
-!! AUTHOR
-!! ------
-!! V. Masson * Meteo France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 13/10/99
-!! V. Masson 06/11/02 optimization and add Businger fonction for TKE
-!! V. Masson 01/01/03 use PAULSON_PSIM function
-!-----------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-!
-!
-INTERFACE BUSINGER_PHIM
- MODULE PROCEDURE BUSINGER_PHIM_0D
- MODULE PROCEDURE BUSINGER_PHIM_1D
- MODULE PROCEDURE BUSINGER_PHIM_2D
- MODULE PROCEDURE BUSINGER_PHIM_3D
-END INTERFACE
-INTERFACE BUSINGER_PHIH
- MODULE PROCEDURE BUSINGER_PHIH_0D
- MODULE PROCEDURE BUSINGER_PHIH_1D
- MODULE PROCEDURE BUSINGER_PHIH_2D
- MODULE PROCEDURE BUSINGER_PHIH_3D
-END INTERFACE
-INTERFACE BUSINGER_PHIE
- MODULE PROCEDURE BUSINGER_PHIE_3D
-END INTERFACE
-INTERFACE PAULSON_PSIM
- MODULE PROCEDURE PAULSON_PSIM_0D
- MODULE PROCEDURE PAULSON_PSIM_1D
- MODULE PROCEDURE PAULSON_PSIM_2D
-END INTERFACE
-INTERFACE LMO
- MODULE PROCEDURE LMO_0D
- MODULE PROCEDURE LMO_1D
- MODULE PROCEDURE LMO_2D
-END INTERFACE
-INTERFACE USTAR
- MODULE PROCEDURE USTAR_0D
- MODULE PROCEDURE USTAR_1D
- MODULE PROCEDURE USTAR_2D
-END INTERFACE
-!
-!-------------------------------------------------------------------------------
-CONTAINS
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIM_3D(PZ_O_LMO)
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1), &
- SIZE(PZ_O_LMO,2),SIZE(PZ_O_LMO,3)) :: BUSINGER_PHIM_3D
-!
- WHERE ( PZ_O_LMO(:,:,:) < 0. )
- BUSINGER_PHIM_3D(:,:,:) = (1.-15.*PZ_O_LMO)**(-0.25)
- ELSEWHERE
- BUSINGER_PHIM_3D(:,:,:) = 1. + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIM_3D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIM_2D(PZ_O_LMO)
- REAL, DIMENSION(:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1),SIZE(PZ_O_LMO,2)) :: BUSINGER_PHIM_2D
-!
- WHERE ( PZ_O_LMO(:,:) < 0. )
- BUSINGER_PHIM_2D(:,:) = (1.-15.*PZ_O_LMO)**(-0.25)
- ELSEWHERE
- BUSINGER_PHIM_2D(:,:) = 1. + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIM_2D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIM_1D(PZ_O_LMO)
- REAL, DIMENSION(:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO)) :: BUSINGER_PHIM_1D
-!
- WHERE ( PZ_O_LMO(:) < 0. )
- BUSINGER_PHIM_1D(:) = (1.-15.*PZ_O_LMO)**(-0.25)
- ELSEWHERE
- BUSINGER_PHIM_1D(:) = 1. + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIM_1D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIM_0D(PZ_O_LMO)
- REAL, INTENT(IN) :: PZ_O_LMO
- REAL :: BUSINGER_PHIM_0D
-!
- IF ( PZ_O_LMO < 0. ) THEN
- BUSINGER_PHIM_0D = (1.-15.*PZ_O_LMO)**(-0.25)
- ELSE
- BUSINGER_PHIM_0D = 1. + 4.7 * PZ_O_LMO
- END IF
-END FUNCTION BUSINGER_PHIM_0D
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIH_3D(PZ_O_LMO)
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1), &
- SIZE(PZ_O_LMO,2),SIZE(PZ_O_LMO,3)) :: BUSINGER_PHIH_3D
-!
- WHERE ( PZ_O_LMO(:,:,:) < 0. )
- BUSINGER_PHIH_3D(:,:,:) = 0.74 * (1.-9.*PZ_O_LMO)**(-0.5)
- ELSEWHERE
- BUSINGER_PHIH_3D(:,:,:) = 0.74 + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIH_3D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIH_2D(PZ_O_LMO)
- REAL, DIMENSION(:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1),SIZE(PZ_O_LMO,2)) :: BUSINGER_PHIH_2D
-!
- WHERE ( PZ_O_LMO(:,:) < 0. )
- BUSINGER_PHIH_2D(:,:) = 0.74 * (1.-9.*PZ_O_LMO)**(-0.5)
- ELSEWHERE
- BUSINGER_PHIH_2D(:,:) = 0.74 + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIH_2D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIH_1D(PZ_O_LMO)
- REAL, DIMENSION(:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO)) :: BUSINGER_PHIH_1D
-!
- WHERE ( PZ_O_LMO(:) < 0. )
- BUSINGER_PHIH_1D(:) = 0.74 * (1.-9.*PZ_O_LMO)**(-0.5)
- ELSEWHERE
- BUSINGER_PHIH_1D(:) = 0.74 + 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION BUSINGER_PHIH_1D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIH_0D(PZ_O_LMO)
- REAL, INTENT(IN) :: PZ_O_LMO
- REAL :: BUSINGER_PHIH_0D
-!
- IF ( PZ_O_LMO < 0. ) THEN
- BUSINGER_PHIH_0D = 0.74 * (1.-9.*PZ_O_LMO)**(-0.5)
- ELSE
- BUSINGER_PHIH_0D = 0.74 + 4.7 * PZ_O_LMO
- END IF
-END FUNCTION BUSINGER_PHIH_0D
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION BUSINGER_PHIE_3D(PZ_O_LMO)
- USE MODD_CTURB
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1), &
- SIZE(PZ_O_LMO,2),SIZE(PZ_O_LMO,3)) :: BUSINGER_PHIE_3D
-!
- WHERE ( PZ_O_LMO(:,:,:) < 0. )
- BUSINGER_PHIE_3D(:,:,:) = (1.+(-PZ_O_LMO)**(2./3.)/XALPSBL) &
- * (1.-15.*PZ_O_LMO)**(0.5)
- ELSEWHERE
- BUSINGER_PHIE_3D(:,:,:) = 1./(1. + 4.7 * PZ_O_LMO)**2
- END WHERE
-END FUNCTION BUSINGER_PHIE_3D
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION PAULSON_PSIM_2D(PZ_O_LMO)
- USE MODD_CST
- REAL, DIMENSION(:,:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1),SIZE(PZ_O_LMO,2)) :: PAULSON_PSIM_2D
-!
- REAL, DIMENSION(SIZE(PZ_O_LMO,1),SIZE(PZ_O_LMO,2)) :: ZX
-
- ZX=1.
- WHERE ( PZ_O_LMO(:,:) < 0. )
- ZX=(1.-15.*PZ_O_LMO)**(0.25)
- PAULSON_PSIM_2D(:,:) = LOG( (1.+ZX**2)*(1+ZX)**2/8. ) - 2.*ATAN(ZX) + XPI/2.
- ELSEWHERE
- PAULSON_PSIM_2D(:,:) = - 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION PAULSON_PSIM_2D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION PAULSON_PSIM_1D(PZ_O_LMO)
- USE MODD_CST
- REAL, DIMENSION(:), INTENT(IN) :: PZ_O_LMO
- REAL, DIMENSION(SIZE(PZ_O_LMO,1)) :: PAULSON_PSIM_1D
-!
- REAL, DIMENSION(SIZE(PZ_O_LMO,1)) :: ZX
-
- ZX=1.
- WHERE ( PZ_O_LMO(:) < 0. )
- ZX=(1.-15.*PZ_O_LMO)**(0.25)
- PAULSON_PSIM_1D(:) = LOG( (1.+ZX**2)*(1+ZX)**2/8. ) - 2.*ATAN(ZX) + XPI/2.
- ELSEWHERE
- PAULSON_PSIM_1D(:) = - 4.7 * PZ_O_LMO
- END WHERE
-END FUNCTION PAULSON_PSIM_1D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION PAULSON_PSIM_0D(PZ_O_LMO)
- USE MODD_CST
- REAL, INTENT(IN) :: PZ_O_LMO
- REAL :: PAULSON_PSIM_0D
-!
- REAL :: ZX
-
- ZX=1.
- IF ( PZ_O_LMO < 0. ) THEN
- ZX=(1.-15.*PZ_O_LMO)**(0.25)
- PAULSON_PSIM_0D = LOG( (1.+ZX**2)*(1+ZX)**2/8. ) - 2.*ATAN(ZX) + XPI/2.
- ELSE
- PAULSON_PSIM_0D = - 4.7 * PZ_O_LMO
- END IF
-END FUNCTION PAULSON_PSIM_0D
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION LMO_2D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, DIMENSION(:,:), INTENT(IN) :: PUSTAR
- REAL, DIMENSION(:,:), INTENT(IN) :: PTHETA
- REAL, DIMENSION(:,:), INTENT(IN) :: PRV
- REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH
- REAL, DIMENSION(:,:), INTENT(IN) :: PSFRV
- REAL, DIMENSION(SIZE(PUSTAR,1),SIZE(PUSTAR,2)) :: LMO_2D
-!
- REAL, DIMENSION(SIZE(PUSTAR,1),SIZE(PUSTAR,2)) :: ZTHETAV
- REAL, DIMENSION(SIZE(PUSTAR,1),SIZE(PUSTAR,2)) :: ZQ0
- REAL :: ZEPS
-!
-!
- ZEPS=(XRV-XRD)/XRD
- ZTHETAV(:,:) = PTHETA(:,:) * ( 1. +ZEPS * PRV(:,:))
- ZQ0 (:,:) = PSFTH(:,:) + ZTHETAV(:,:) * ZEPS * PSFRV(:,:)
-!
- LMO_2D(:,:) = XUNDEF
- WHERE ( ZQ0(:,:) /=0. ) &
- LMO_2D(:,:) = - MAX(PUSTAR(:,:),1.E-6)**3 &
- / ( XKARMAN * XG / ZTHETAV(:,:) *ZQ0(:,:) )
-
-END FUNCTION LMO_2D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION LMO_1D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, DIMENSION(:), INTENT(IN) :: PUSTAR
- REAL, DIMENSION(:), INTENT(IN) :: PTHETA
- REAL, DIMENSION(:), INTENT(IN) :: PRV
- REAL, DIMENSION(:), INTENT(IN) :: PSFTH
- REAL, DIMENSION(:), INTENT(IN) :: PSFRV
- REAL, DIMENSION(SIZE(PUSTAR)) :: LMO_1D
-!
- REAL, DIMENSION(SIZE(PUSTAR)) :: ZTHETAV
- REAL :: ZEPS
-!
-!
- ZEPS=(XRV-XRD)/XRD
-!
- ZTHETAV(:) = PTHETA(:) * ( 1. +ZEPS * PRV(:))
-!
- LMO_1D(:) = XUNDEF
- WHERE ( PSFTH(:)/ZTHETAV(:)+ZEPS*PSFRV(:)/=0. ) &
- LMO_1D(:) = - MAX(PUSTAR(:),1.E-6)**3 &
- / ( XKARMAN * ( XG / ZTHETAV(:) * PSFTH(:) &
- + XG * ZEPS * PSFRV(:) ) )
-END FUNCTION LMO_1D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION LMO_0D(PUSTAR,PTHETA,PRV,PSFTH,PSFRV)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, INTENT(IN) :: PUSTAR
- REAL, INTENT(IN) :: PTHETA
- REAL, INTENT(IN) :: PRV
- REAL, INTENT(IN) :: PSFTH
- REAL, INTENT(IN) :: PSFRV
- REAL :: LMO_0D
-!
- REAL :: ZTHETAV
- REAL :: ZEPS
-!
-!
- ZEPS=(XRV-XRD)/XRD
-!
-!
- ZTHETAV = PTHETA * ( 1. +ZEPS * PRV)
-!
- LMO_0D = XUNDEF
- IF ( PSFTH/ZTHETAV+ZEPS*PSFRV/=0. ) &
- LMO_0D = - MAX(PUSTAR,1.E-6)**3 &
- / ( XKARMAN * ( XG / ZTHETAV * PSFTH &
- + XG * ZEPS * PSFRV ) )
-END FUNCTION LMO_0D
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-FUNCTION USTAR_2D(PU,PV,PZ,PZ0,PLMO)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, DIMENSION(:,:), INTENT(IN) :: PU
- REAL, DIMENSION(:,:), INTENT(IN) :: PV
- REAL, DIMENSION(:,:), INTENT(IN) :: PZ
- REAL, DIMENSION(:,:), INTENT(IN) :: PZ0
- REAL, DIMENSION(:,:), INTENT(IN) :: PLMO
- REAL, DIMENSION(SIZE(PU,1),SIZE(PU,2)) :: USTAR_2D
-
- REAL, DIMENSION(SIZE(PU,1),SIZE(PU,2)) :: ZZ_O_LMO
- REAL, DIMENSION(SIZE(PU,1),SIZE(PU,2)) :: ZZ0_O_LMO
-!
-!* purely unstable case
- USTAR_2D(:,:) = 0.
- ZZ_O_LMO(:,:) = XUNDEF
- ZZ0_O_LMO(:,:) = XUNDEF
-!
-!* general case
- WHERE(ABS(PLMO) > 1.E-20 .AND. PLMO/=XUNDEF)
- ZZ_O_LMO = PZ(:,:) / PLMO(:,:)
- ZZ0_O_LMO = PZ0(:,:) / PLMO(:,:)
- USTAR_2D(:,:) = SQRT( PU(:,:)**2+PV(:,:)**2 ) &
- * XKARMAN / ( LOG(PZ(:,:)/PZ0(:,:)) &
- - PAULSON_PSIM(ZZ_O_LMO(:,:)) &
- + PAULSON_PSIM(ZZ0_O_LMO(:,:)) )
- END WHERE
-!
-!* purely neutral case
- WHERE(PLMO==XUNDEF)
- ZZ_O_LMO = 0.
- USTAR_2D(:,:) = SQRT( PU(:,:)**2+PV(:,:)**2 ) &
- * XKARMAN / LOG(PZ(:,:)/PZ0(:,:))
- END WHERE
-!
-END FUNCTION USTAR_2D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION USTAR_1D(PU,PV,PZ,PZ0,PLMO)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, DIMENSION(:), INTENT(IN) :: PU
- REAL, DIMENSION(:), INTENT(IN) :: PV
- REAL, DIMENSION(:), INTENT(IN) :: PZ
- REAL, DIMENSION(:), INTENT(IN) :: PZ0
- REAL, DIMENSION(:), INTENT(IN) :: PLMO
- REAL, DIMENSION(SIZE(PU)) :: USTAR_1D
-
- REAL, DIMENSION(SIZE(PU)) :: ZZ_O_LMO
- REAL, DIMENSION(SIZE(PU)) :: ZZ0_O_LMO
-!
-!* purely unstable case
- USTAR_1D(:) = 0.
- ZZ_O_LMO(:) = XUNDEF
- ZZ0_O_LMO(:) = XUNDEF
-!
-!* general case
- WHERE(ABS(PLMO) > 1.E-20 .AND. PLMO/=XUNDEF)
- ZZ_O_LMO = PZ(:) / PLMO(:)
- ZZ0_O_LMO = PZ0(:) / PLMO(:)
- USTAR_1D(:) = SQRT( PU(:)**2+PV(:)**2 ) &
- * XKARMAN / ( LOG(PZ(:)/PZ0(:)) &
- - PAULSON_PSIM(ZZ_O_LMO(:)) &
- + PAULSON_PSIM(ZZ0_O_LMO(:)) )
- END WHERE
-!
-!* purely neutral case
- WHERE(PLMO==XUNDEF)
- ZZ_O_LMO = 0.
- USTAR_1D(:) = SQRT( PU(:)**2+PV(:)**2 ) &
- * XKARMAN / LOG(PZ(:)/PZ0(:))
- END WHERE
-!
-END FUNCTION USTAR_1D
-!
-!-------------------------------------------------------------------------------
-!
-FUNCTION USTAR_0D(PU,PV,PZ,PZ0,PLMO)
- USE MODD_CST
- USE MODD_PARAMETERS
- REAL, INTENT(IN) :: PU
- REAL, INTENT(IN) :: PV
- REAL, INTENT(IN) :: PZ
- REAL, INTENT(IN) :: PZ0
- REAL, INTENT(IN) :: PLMO
- REAL :: USTAR_0D
-!
-!* purely unstable case
- USTAR_0D = 0.
-!
-!* general case
- IF ( ABS(PLMO) >= 1.E-20 .AND. PLMO/=XUNDEF) &
- USTAR_0D = SQRT( PU**2+PV**2 ) &
- * XKARMAN / ( LOG(PZ/PZ0) &
- - PAULSON_PSIM(PZ/PLMO) &
- + PAULSON_PSIM(PZ0/PLMO))
-!
-!* purely neutral case
- IF (PLMO==XUNDEF) &
- USTAR_0D = SQRT( PU**2+PV**2 ) &
- * XKARMAN / LOG(PZ/PZ0)
-
-END FUNCTION USTAR_0D
-!
-!-------------------------------------------------------------------------------
-!
-END MODULE MODE_SBL