diff --git a/src/arome/aux/mode_sources_neg_correct.F90 b/src/arome/aux/mode_sources_neg_correct.F90
new file mode 100644
index 0000000000000000000000000000000000000000..49bf5f75ddcb4a854bb3ddde0bd2c3c37c47be4c
--- /dev/null
+++ b/src/arome/aux/mode_sources_neg_correct.F90
@@ -0,0 +1,22 @@
+MODULE MODE_SOURCES_NEG_CORRECT
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE SOURCES_NEG_CORRECT_PHY(D, KSV, HCLOUD, HBUDNAME, KRR, PTSTEP, PPABST, &
+ &PTHT, PRT, PRTHS, PRRS, PRSVS, PRHODJ)
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+IMPLICIT NONE
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KSV ! Number of SV variables
+CHARACTER(LEN=*), INTENT(IN) :: HCLOUD ! Kind of cloud parameterization
+CHARACTER(LEN=*), INTENT(IN) :: HBUDNAME ! Budget name
+INTEGER, INTENT(IN) :: KRR ! Number of moist variables
+REAL, INTENT(IN) :: PTSTEP ! Timestep
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PPABST ! Absolute pressure at time t
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), INTENT(IN) :: PRT ! Moist variables at time t
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRTHS ! Source terms
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), INTENT(INOUT) :: PRRS ! Source terms
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KSV), INTENT(INOUT) :: PRSVS ! Source terms
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN), OPTIONAL :: PRHODJ ! Dry density * jacobian
+END SUBROUTINE SOURCES_NEG_CORRECT_PHY
+END MODULE MODE_SOURCES_NEG_CORRECT
diff --git a/src/common/aux/ibm_mixinglength.F90 b/src/common/aux/ibm_mixinglength.F90
deleted file mode 100644
index 766627888e3d49a472bd2eaa7a0b49bb58667abe..0000000000000000000000000000000000000000
--- a/src/common/aux/ibm_mixinglength.F90
+++ /dev/null
@@ -1,35 +0,0 @@
-!MNH_LIC Copyright 2019-2021 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 MODI_IBM_MIXINGLENGTH
- ! ############################
- !
- INTERFACE
- !
- SUBROUTINE IBM_MIXINGLENGTH(PLM,PLEPS,PMU,PHI,PTKE)
- !
- REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM
- REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS
- REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMU
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PHI
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
- !
- END SUBROUTINE IBM_MIXINGLENGTH
- !
- END INTERFACE
-END MODULE MODI_IBM_MIXINGLENGTH
- !
- SUBROUTINE IBM_MIXINGLENGTH(PLM,PLEPS,PMU,PHI,PTKE)
- !
- REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM
- REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS
- REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMU
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PHI
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
- !
- END SUBROUTINE IBM_MIXINGLENGTH
- !
diff --git a/src/common/aux/mode_sources_neg_correct.F90 b/src/common/aux/mode_sources_neg_correct.F90
deleted file mode 100644
index 1b49a6e2b22a64880bea7e75e5c32001a106ecda..0000000000000000000000000000000000000000
--- a/src/common/aux/mode_sources_neg_correct.F90
+++ /dev/null
@@ -1,19 +0,0 @@
-MODULE MODE_SOURCES_NEG_CORRECT
-IMPLICIT NONE
-CONTAINS
-SUBROUTINE SOURCES_NEG_CORRECT(HCLOUD, HBUDNAME, KRR, PTSTEP, PPABST, &
- &PTHT, PRT, PRTHS, PRRS, PRSVS, PRHODJ)
-IMPLICIT NONE
-CHARACTER(LEN=*), INTENT(IN) :: HCLOUD ! Kind of cloud parameterization
-CHARACTER(LEN=*), INTENT(IN) :: HBUDNAME ! Budget name
-INTEGER, INTENT(IN) :: KRR ! Number of moist variables
-REAL, INTENT(IN) :: PTSTEP ! Timestep
-REAL, DIMENSION(:, :, :), INTENT(IN) :: PPABST ! Absolute pressure at time t
-REAL, DIMENSION(:, :, :), INTENT(IN) :: PTHT ! Theta at time t
-REAL, DIMENSION(:, :, :, :), INTENT(IN) :: PRT ! Moist variables at time t
-REAL, DIMENSION(:, :, :), INTENT(INOUT) :: PRTHS ! Source terms
-REAL, DIMENSION(:, :, :, :), INTENT(INOUT) :: PRRS ! Source terms
-REAL, DIMENSION(:, :, :, :), INTENT(INOUT) :: PRSVS ! Source terms
-REAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: PRHODJ ! Dry density * jacobian
-END SUBROUTINE SOURCES_NEG_CORRECT
-END MODULE MODE_SOURCES_NEG_CORRECT
diff --git a/src/common/turb/mode_ibm_mixinglength.F90 b/src/common/turb/mode_ibm_mixinglength.F90
new file mode 100644
index 0000000000000000000000000000000000000000..7f74c571a60118414f48f509d1f3b7f95e3d5e26
--- /dev/null
+++ b/src/common/turb/mode_ibm_mixinglength.F90
@@ -0,0 +1,146 @@
+!MNH_LIC Copyright 2019-2021 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_IBM_MIXINGLENGTH
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE IBM_MIXINGLENGTH(D,PLM,PLEPS,PMU,PHI,PTKE)
+ ! ###################################################
+ !
+ !**** *IBM_MIXINGLENGTH* - Alteration of the mixing lenght (IBM)
+ !
+ ! PURPOSE
+ ! -------
+ ! The limitation is corrected for the immersed bonudary method:
+ ! => using the level set phi
+ ! => LM < k(-phi)
+ !
+ ! METHOD
+ ! ------
+ !
+ ! INDEX
+ ! -----
+ !
+ ! IMPLICIT ARGUMENTS
+ ! ------------------
+ !
+ ! REFERENCE
+ ! ---------
+ !
+ ! AUTHOR
+ ! ------
+ !
+ ! Franck Auguste * CERFACS(AE) *
+ !
+ ! MODIFICATIONS
+ ! -------------
+ ! Original 01/01/2019
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !**** 0. DECLARATIONS
+ ! ------------------
+ !
+ ! module
+ USE MODE_ll
+ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+ !
+ ! declaration
+ USE MODD_FIELD_n
+ USE MODD_PARAMETERS
+ USE MODD_IBM_PARAM_n
+ USE MODD_REF_n, ONLY: XRHODJ,XRHODREF
+ USE MODD_CTURB
+ USE MODD_CST
+ USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ
+ !
+ ! interface
+ !
+ IMPLICIT NONE
+ !
+ !------------------------------------------------------------------------------
+ !
+ ! 0.1 Declaration of arguments
+ TYPE(DIMPHYEX_t), INTENT(IN) :: D
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PLM
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PLEPS
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMU
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PHI
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
+ !
+ !------------------------------------------------------------------------------
+ !
+ ! 0.2 Declaration of local variables
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZALPHA,ZBETA
+ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZLM,ZMU,ZLN
+ INTEGER :: IKU,IKB,IKE,IIB,IIE,IJB,IJE
+ REAL :: ZKARMAN
+ !
+ !-------------------------------------------------------------------------------
+ !
+ IKU = D%NKT
+ IKE = D%NKE
+ IKB = D%NKB
+ IIB = D%NIB
+ IIE = D%NIE
+ IJB = D%NJB
+ IJE = D%NJE
+ !
+ ! Turbulent velocity
+ !
+ ZMU(:,:,:)= 2.*XIBM_CNU**0.25*(PTKE(:,:,:))**(1./2.) !2 correspond to KTKE
+ ZKARMAN = XKARMAN*XCED/XIBM_CNU**0.75
+ !
+ ! Mesh scale
+ !
+ ZLN(:,:,:) = (XRHODJ(:,:,:)/XRHODREF(:,:,:))**(1./3.)
+ ZLM(:,:,:) = PLM(:,:,:)
+ !
+ ! limit domain
+ !
+ ZBETA(:,:,:)= XZZ (:,:,:)
+ ZBETA(:,:,IKB:IKE) = 0.5*(XZZ(:,:,IKB+1:IKE+1)+XZZ(:,:,IKB:IKE))
+ ZBETA(:,:,IKB-1) = -0.5*(XZZ(:,:,IKB+1)+XZZ(:,:,IKB))
+ ZBETA(:,:,IKE+1) = ZBETA(:,:,IKE)
+ ZLM(:,:,:) = MIN(ZLM(:,:,:),+ZKARMAN*ZBETA(:,:,:))
+ !
+ ! limit immersed wall
+ !
+ ZLM(:,:,:) = MIN(ZLM(:,:,:),-ZKARMAN*PHI(:,:,:))
+ !
+ ! limit physical scale
+ ZALPHA(:,:,:) = MIN(9.8*XIBM_RUG,0.5*ZKARMAN*ZLN(:,:,:))
+ ZLM(:,:,:) = MAX(ZALPHA(:,:,:),ZLM(:,:,:))
+ !
+ ! Boundary condition
+ ZMU(:,:,IKB-1)= ZMU(:,:,IKB)
+ ZLM(:,:,IKB-1)= ZLM(:,:,IKB)
+ ZMU(:,:,IKE+1)= ZMU(:,:,IKE)
+ ZLM(:,:,IKE+1)= ZLM(:,:,IKE)
+ IF (LEAST_ll()) THEN
+ ZMU(IIE+1,:,:)= ZMU(IIE,:,:)
+ ZLM(IIE+1,:,:)= ZLM(IIE,:,:)
+ ENDIF
+ IF (LWEST_ll()) THEN
+ ZMU(IIB-1,:,:)= ZMU(IIB,:,:)
+ ZLM(IIB-1,:,:)= ZLM(IIB,:,:)
+ ENDIF
+ IF (LNORTH_ll()) THEN
+ ZMU(:,IJE+1,:)= ZMU(:,IJE,:)
+ ZLM(:,IJE+1,:)= ZLM(:,IJE,:)
+ ENDIF
+ IF (LSOUTH_ll()) THEN
+ ZMU(:,IJB-1,:)= ZMU(:,IJB,:)
+ ZLM(:,IJB-1,:)= ZLM(:,IJB,:)
+ ENDIF
+ !
+ !Communication
+ PLM(:,:,:) = ZLM(:,:,:)
+ PLEPS(:,:,:) = PLM(:,:,:)
+ PMU(:,:,:) = ZMU(:,:,:)
+ !
+END SUBROUTINE IBM_MIXINGLENGTH
+END MODULE MODE_IBM_MIXINGLENGTH
diff --git a/src/common/turb/mode_rotate_wind.F90 b/src/common/turb/mode_rotate_wind.F90
index e91e2a20684934c1d88905d62c3da628b10258ea..74c4932ffa5656f1d643b37efd7fbbf41dfee7a0 100644
--- a/src/common/turb/mode_rotate_wind.F90
+++ b/src/common/turb/mode_rotate_wind.F90
@@ -8,7 +8,7 @@
IMPLICIT NONE
CONTAINS
! ###########################################################
- SUBROUTINE ROTATE_WIND(PU,PV,PW, &
+ SUBROUTINE ROTATE_WIND(D,PU,PV,PW, &
PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
PCOSSLOPE,PSINSLOPE, &
PDXX,PDYY,PDZZ, &
@@ -74,34 +74,36 @@ CONTAINS
!* 0. DECLARATIONS
! ------------
USE MODD_PARAMETERS, ONLY: JPVEXT
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
!
IMPLICIT NONE
!
!
!* 0.1 declarations of arguments
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PU,PV,PW ! cartesian components
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PU,PV,PW ! cartesian components
! of the wind
-REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
! Director Cosinus along x, y and z directions at surface w-point
-REAL, DIMENSION(:,:), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
! between i and the slope vector
-REAL, DIMENSION(:,:), INTENT(IN) :: PSINSLOPE ! sinus of the angle
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(IN) :: PSINSLOPE ! sinus of the angle
! between i and the slope vector
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDXX, PDYY, PDZZ
! Metric coefficients
-REAL, DIMENSION(:,:), INTENT(OUT) :: PUSLOPE ! wind component along
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PUSLOPE ! wind component along
! the maximum slope direction
-REAL, DIMENSION(:,:), INTENT(OUT) :: PVSLOPE ! wind component along
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PVSLOPE ! wind component along
! the direction normal to the maximum slope one
!
!-------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
!
-INTEGER, DIMENSION(SIZE(PDIRCOSXW,1),SIZE(PDIRCOSXW,2)) :: ILOC,JLOC
+INTEGER, DIMENSION(D%NIT,D%NJT) :: ILOC,JLOC
! shift index to find the 4 nearest points in x and y directions
-REAL, DIMENSION(SIZE(PDIRCOSXW,1),SIZE(PDIRCOSXW,2)) :: ZCOEFF,ZCOEFM, &
+REAL, DIMENSION(D%NIT,D%NJT) :: ZCOEFF,ZCOEFM, &
! interpolation weigths for flux and mass locations
ZUINT,ZVINT,ZWINT, &
! intermediate values of the cartesian components after x interp.
@@ -121,8 +123,8 @@ INTEGER :: JI,JJ
!* 1. PRELIMINARIES
! -------------
!
-PUSLOPE=0.
-PVSLOPE=0.
+PUSLOPE(:,:)=0.
+PVSLOPE(:,:)=0.
!
IIB = 2
IJB = 2
@@ -201,4 +203,76 @@ END DO
!----------------------------------------------------------------------------
!
END SUBROUTINE ROTATE_WIND
+!
+! ##############################################
+ SUBROUTINE UPDATE_ROTATE_WIND(D,PUSLOPE,PVSLOPE,HLBCX,HLBCY)
+! ##############################################
+!!
+!!**** *UPDATE_ROTATE_WIND* routine to set rotate wind values at the border
+!
+!! AUTHOR
+!! ------
+!!
+!! P Jabouille *CNRM METEO-FRANCE
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/06/99
+!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODE_ll, ONLY: ADD2DFIELD_ll, UPDATE_HALO_ll, CLEANLIST_ll, &
+ LWEST_ll, LEAST_ll, LSOUTH_ll, LNORTH_ll
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
+REAL, DIMENSION(D%NIT,D%NJT), INTENT(INOUT) :: PUSLOPE,PVSLOPE
+!
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+INTEGER :: IINFO_ll ! return code of parallel routine
+!
+! tangential surface fluxes in the axes following the orography
+!
+!* 1 PROLOGUE
+!
+NULLIFY(TZFIELDS_ll)
+!
+! 2 Update halo if necessary
+!
+!!$IF (NHALO == 1) THEN
+ CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' )
+ CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+!!$ENDIF
+!
+! 3 Boundary conditions for non cyclic case
+!
+IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
+ PUSLOPE(D%NIB-1,:)=PUSLOPE(D%NIB,:)
+ PVSLOPE(D%NIB-1,:)=PVSLOPE(D%NIB,:)
+END IF
+IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
+ PUSLOPE(D%NIE+1,:)=PUSLOPE(D%NIE,:)
+ PVSLOPE(D%NIE+1,:)=PVSLOPE(D%NIE,:)
+END IF
+IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
+ PUSLOPE(:,D%NJB-1)=PUSLOPE(:,D%NJB)
+ PVSLOPE(:,D%NJB-1)=PVSLOPE(:,D%NJB)
+END IF
+IF( HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
+ PUSLOPE(:,D%NJE+1)=PUSLOPE(:,D%NJE)
+ PVSLOPE(:,D%NJE+1)=PVSLOPE(:,D%NJE)
+END IF
+!
+END SUBROUTINE UPDATE_ROTATE_WIND
END MODULE MODE_ROTATE_WIND
diff --git a/src/common/turb/mode_sbl_phy.F90 b/src/common/turb/mode_sbl_phy.F90
index a68007cbf958791bea23ff491b034dd0e2e2b205..a45b6485a35706a1397703c31ae0b7d179c21d96 100644
--- a/src/common/turb/mode_sbl_phy.F90
+++ b/src/common/turb/mode_sbl_phy.F90
@@ -58,16 +58,20 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZ_O_LMO
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: BUSINGERPHIM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER :: JIJ,JK
+INTEGER :: JIJ,JK,IIJB,IIJE
!
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM',0,ZHOOK_HANDLE)
-!$mnh_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
-WHERE ( PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT) < 0. )
- BUSINGERPHIM(1:D%NIEC*D%NJEC,1:D%NKT) = (1.-15.*PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT))**(-0.25)
+!
+IIJE=D%NIJE
+IIJB=D%NIJB
+!
+!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+WHERE ( PZ_O_LMO(IIJB:IIJE,1:D%NKT) < 0. )
+ BUSINGERPHIM(IIJB:IIJE,1:D%NKT) = (1.-15.*PZ_O_LMO(IIJB:IIJE,1:D%NKT))**(-0.25)
ELSEWHERE
- BUSINGERPHIM(1:D%NIEC*D%NJEC,1:D%NKT) = 1. + 4.7 * PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT)
+ BUSINGERPHIM(IIJB:IIJE,1:D%NKT) = 1. + 4.7 * PZ_O_LMO(IIJB:IIJE,1:D%NKT)
END WHERE
-!$mnh_end_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
+!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIM',1,ZHOOK_HANDLE)
END SUBROUTINE BUSINGER_PHIM
!
@@ -84,16 +88,20 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZ_O_LMO
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: BUSINGERPHIH
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER :: JIJ,JK
+INTEGER :: JIJ,JK,IIJB,IIJE
!
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH',0,ZHOOK_HANDLE)
-!$mnh_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
-WHERE ( PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT) < 0. )
- BUSINGERPHIH(1:D%NIEC*D%NJEC,1:D%NKT) = 0.74 * (1.-9.*PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT))**(-0.5)
+!
+IIJE=D%NIJE
+IIJB=D%NIJB
+!
+!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+WHERE ( PZ_O_LMO(IIJB:IIJE,1:D%NKT) < 0. )
+ BUSINGERPHIH(IIJB:IIJE,1:D%NKT) = 0.74 * (1.-9.*PZ_O_LMO(IIJB:IIJE,1:D%NKT))**(-0.5)
ELSEWHERE
- BUSINGERPHIH(1:D%NIEC*D%NJEC,1:D%NKT) = 0.74 + 4.7 * PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT)
+ BUSINGERPHIH(IIJB:IIJE,1:D%NKT) = 0.74 + 4.7 * PZ_O_LMO(IIJB:IIJE,1:D%NKT)
END WHERE
-!$mnh_end_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
+!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIH',1,ZHOOK_HANDLE)
END SUBROUTINE BUSINGER_PHIH
!
@@ -111,17 +119,62 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZ_O_LMO
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT) :: BUSINGERPHIE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER :: JIJ,JK
+INTEGER :: JIJ,JK,IIJB,IIJE
!
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIE',0,ZHOOK_HANDLE)
-!$mnh_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
-WHERE ( PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT) < 0. )
- BUSINGERPHIE(1:D%NIEC*D%NJEC,1:D%NKT)=(1.+(-PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT))**(2./3.)/CSTURB%XALPSBL)&
- * (1.-15.*PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT))**(0.5)
+!
+IIJE=D%NIJE
+IIJB=D%NIJB
+!
+!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+WHERE ( PZ_O_LMO(IIJB:IIJE,1:D%NKT) < 0. )
+ BUSINGERPHIE(IIJB:IIJE,1:D%NKT)=(1.+(-PZ_O_LMO(IIJB:IIJE,1:D%NKT))**(2./3.)/CSTURB%XALPSBL)&
+ * (1.-15.*PZ_O_LMO(IIJB:IIJE,1:D%NKT))**(0.5)
ELSEWHERE
- BUSINGERPHIE(1:D%NIEC*D%NJEC,1:D%NKT) = 1./(1. + 4.7 * PZ_O_LMO(1:D%NIEC*D%NJEC,1:D%NKT))**2
+ BUSINGERPHIE(IIJB:IIJE,1:D%NKT) = 1./(1. + 4.7 * PZ_O_LMO(IIJB:IIJE,1:D%NKT))**2
END WHERE
-!$mnh_end_expand_where(JIJ=1:D%NIEC*D%NJEC,JK=1:D%NKT)
+!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
IF (LHOOK) CALL DR_HOOK('MODE_SBL:BUSINGER_PHIE',1,ZHOOK_HANDLE)
END SUBROUTINE BUSINGER_PHIE
+!
+SUBROUTINE LMO(D,CST,PUSTAR,PTHETA,PRV,PSFTH,PSFRV,PLMO)
+ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+ USE MODD_CST, ONLY: CST_t
+ USE MODD_PARAMETERS, ONLY: JPVEXT_TURB,XUNDEF
+ !
+ TYPE(DIMPHYEX_t), INTENT(IN) :: D
+ TYPE(CST_t), INTENT(IN) :: CST
+ REAL, DIMENSION(D%NIJT), INTENT(IN) :: PUSTAR
+ REAL, DIMENSION(D%NIJT), INTENT(IN) :: PTHETA
+ REAL, DIMENSION(D%NIJT), INTENT(IN) :: PRV
+ REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFTH
+ REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFRV
+ REAL, DIMENSION(D%NIJT),INTENT(OUT) :: PLMO
+!
+ REAL, DIMENSION(D%NIJT) :: ZTHETAV, ZQ0
+ REAL :: ZEPS
+ INTEGER :: IIJB,IIJE, JIJ
+!
+ REAL(KIND=JPRB) :: ZHOOK_HANDLE
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO',0,ZHOOK_HANDLE)
+!
+ IIJE=D%NIJE
+ IIJB=D%NIJB
+ ZEPS=(CST%XRV-CST%XRD)/CST%XRD
+!
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZTHETAV(IIJB:IIJE) = PTHETA(IIJB:IIJE) * ( 1. +ZEPS * PRV(IIJB:IIJE))
+ ZQ0(IIJB:IIJE) = PSFTH(IIJB:IIJE) + ZTHETAV(IIJB:IIJE) * ZEPS * PSFRV(IIJB:IIJE)
+!
+ PLMO(IIJB:IIJE) = XUNDEF
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ !$mnh_expand_where(JIJ=IIJB:IIJE)
+ WHERE ( ZQ0(IIJB:IIJE)/=0. )
+ PLMO(IIJB:IIJE) = - MAX(PUSTAR(IIJB:IIJE),1.E-6)**3 &
+ / ( CST%XKARMAN * CST%XG / ZTHETAV(IIJB:IIJE) *ZQ0(IIJB:IIJE) )
+ END WHERE
+ !$mnh_end_expand_where(JIJ=IIJB:IIJE)
+ IF (LHOOK) CALL DR_HOOK('MODE_SBL:LMO',1,ZHOOK_HANDLE)
+END SUBROUTINE LMO
+!
END MODULE MODE_SBL_PHY
diff --git a/src/common/turb/mode_tm06.F90 b/src/common/turb/mode_tm06.F90
index 903ea792e01c5f456289a0c283f6e24c76400a41..d6ea7a26875087c1d186639f282a9d10b027d1f8 100644
--- a/src/common/turb/mode_tm06.F90
+++ b/src/common/turb/mode_tm06.F90
@@ -5,7 +5,7 @@
MODULE MODE_TM06
IMPLICIT NONE
CONTAINS
-SUBROUTINE TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,PMWTH,PMTH2)
+SUBROUTINE TM06(D,CST,PTHVREF,PBL_DEPTH,PZZ,PSFTH,PMWTH,PMTH2)
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! #################################################################
@@ -45,10 +45,9 @@ SUBROUTINE TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,PMWTH,PMTH2)
!* 0. DECLARATIONS
! ------------
!
-USE MODD_PARAMETERS, ONLY : XUNDEF
-USE MODD_CST, ONLY : XG
-USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
-
+USE MODD_CST, ONLY: CST_t
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS, ONLY: XUNDEF,JPVEXT_TURB
!
!
IMPLICIT NONE
@@ -56,84 +55,98 @@ IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
-INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
-INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! reference potential temperature
-REAL, DIMENSION(:,:), INTENT(IN) :: PBL_DEPTH ! boundary layer height
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux levels
-REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH ! surface heat flux
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PMWTH ! w'2th'
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PMTH2 ! w'th'2
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+TYPE(CST_t), INTENT(IN) :: CST
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHVREF ! reference potential temperature
+REAL, DIMENSION(D%NIJT), INTENT(IN) :: PBL_DEPTH ! boundary layer height
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZZ ! altitude of flux levels
+REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFTH ! surface heat flux
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT):: PMWTH ! w'2th'
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT):: PMTH2 ! w'th'2
!
!-------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
!
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZZ_O_H ! normalized height z/h (where h=BL height)
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZWSTAR ! normalized convective velocity w*
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZTSTAR ! normalized temperature velocity w*
+REAL, DIMENSION(D%NIJT,D%NKT):: ZZ_O_H ! normalized height z/h (where h=BL height)
+REAL, DIMENSION(D%NIJT) :: ZWSTAR ! normalized convective velocity w*
+REAL, DIMENSION(D%NIJT) :: ZTSTAR ! normalized temperature velocity w*
!
INTEGER :: JK ! loop counter
-INTEGER :: IKT ! vertical size
-INTEGER :: IKTB,IKTE,IKB,IKE ! vertical levels
+INTEGER :: IIJE,IIJB
+INTEGER :: IKTB,IKTE,IKB,IKE,IKT ! vertical levels
!----------------------------------------------------------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TM06',0,ZHOOK_HANDLE)
-IKT=SIZE(PZZ,3)
-IKTB=1+JPVEXT_TURB
-IKTE=IKT-JPVEXT_TURB
-IKB=KKA+JPVEXT_TURB*KKL
-IKE=KKU-JPVEXT_TURB*KKL
-
+IKTB=D%NKTB
+IKTE=D%NKTE
+IKB=D%NKB
+IKT=D%NKT
+IKE=D%NKE
+IIJE=D%NIJE
+IIJB=D%NIJB
!
!
!* w* and T*
!
-WHERE(PSFTH>0.)
- ZWSTAR = ((XG/PTHVREF(:,:,IKB))*PSFTH*PBL_DEPTH)**(1./3.)
- ZTSTAR = PSFTH / ZWSTAR
+!$mnh_expand_where(JIJ=IIJB:IIJE)
+WHERE(PSFTH(IIJB:IIJE)>0.)
+ ZWSTAR(IIJB:IIJE) = ((CST%XG/PTHVREF(IIJB:IIJE,IKB))*PSFTH(IIJB:IIJE)*PBL_DEPTH(IIJB:IIJE))**(1./3.)
+ ZTSTAR(IIJB:IIJE) = PSFTH(IIJB:IIJE) / ZWSTAR(IIJB:IIJE)
ELSEWHERE
- ZWSTAR = 0.
- ZTSTAR = 0.
+ ZWSTAR(IIJB:IIJE) = 0.
+ ZTSTAR(IIJB:IIJE) = 0.
END WHERE
+!$mnh_end_expand_where(JIJ=IIJB:IIJE)
!
!
!* normalized height
!
-ZZ_O_H = XUNDEF
+!$mnh_expand_array(JIJ=IIJB:IIJE)
+ZZ_O_H(IIJB:IIJE,1:D%NKT) = XUNDEF
+!$mnh_end_expand_array(JIJ=IIJB:IIJE)
DO JK=1,IKT
- WHERE (PBL_DEPTH/=XUNDEF)
- ZZ_O_H(:,:,JK) = (PZZ(:,:,JK)-PZZ(:,:,IKB)) / PBL_DEPTH(:,:)
+ !$mnh_expand_where(JIJ=IIJB:IIJE)
+ WHERE (PBL_DEPTH(IIJB:IIJE)/=XUNDEF)
+ ZZ_O_H(IIJB:IIJE,JK) = (PZZ(IIJB:IIJE,JK)-PZZ(IIJB:IIJE,IKB)) / PBL_DEPTH(IIJB:IIJE)
END WHERE
+ !$mnh_end_expand_where(JIJ=IIJB:IIJE)
END DO
!
!* w'th'2
!
-PMTH2 = 0.
-WHERE(ZZ_O_H < 0.95 .AND. ZZ_O_H/=XUNDEF)
- PMTH2(:,:,:) = 4.*(MAX(ZZ_O_H,0.))**0.4*(ZZ_O_H-0.95)**2
+PMTH2(IIJB:IIJE,1:D%NKT) = 0.
+!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+WHERE(ZZ_O_H(IIJB:IIJE,1:D%NKT) < 0.95 .AND. ZZ_O_H(IIJB:IIJE,1:D%NKT)/=XUNDEF)
+ PMTH2(IIJB:IIJE,1:D%NKT) = 4.*(MAX(ZZ_O_H(IIJB:IIJE,1:D%NKT),0.))**0.4*(ZZ_O_H(IIJB:IIJE,1:D%NKT)-0.95)**2
END WHERE
+!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+!$mnh_expand_array(JIJ=IIJB:IIJE)
DO JK=IKTB+1,IKTE-1
- PMTH2(:,:,JK) = PMTH2(:,:,JK) * ZTSTAR(:,:)**2*ZWSTAR(:,:)
+ PMTH2(IIJB:IIJE,JK) = PMTH2(IIJB:IIJE,JK) * ZTSTAR(IIJB:IIJE)**2*ZWSTAR(IIJB:IIJE)
END DO
-PMTH2(:,:,IKE)=PMTH2(:,:,IKE) * ZTSTAR(:,:)**2*ZWSTAR(:,:)
-PMTH2(:,:,KKU)=PMTH2(:,:,KKU) * ZTSTAR(:,:)**2*ZWSTAR(:,:)
-
+PMTH2(IIJB:IIJE,IKE)=PMTH2(IIJB:IIJE,IKE) * ZTSTAR(IIJB:IIJE)**2*ZWSTAR(IIJB:IIJE)
+PMTH2(IIJB:IIJE,D%NKU)=PMTH2(IIJB:IIJE,D%NKU) * ZTSTAR(IIJB:IIJE)**2*ZWSTAR(IIJB:IIJE)
+!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
!
!* w'2th'
!
-PMWTH = 0.
-WHERE(ZZ_O_H <0.9 .AND. ZZ_O_H/=XUNDEF)
- PMWTH(:,:,:) = MAX(-7.9*(ABS(ZZ_O_H-0.35))**2.9 * (ABS(ZZ_O_H-1.))**0.58 + 0.37, 0.)
+PMWTH(IIJB:IIJE,1:D%NKT) = 0.
+!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+WHERE(ZZ_O_H(IIJB:IIJE,1:D%NKT) <0.9 .AND. ZZ_O_H(IIJB:IIJE,1:D%NKT)/=XUNDEF)
+ PMWTH(IIJB:IIJE,1:D%NKT) = MAX(-7.9*(ABS(ZZ_O_H(IIJB:IIJE,1:D%NKT)-0.35))**2.9 &
+ * (ABS(ZZ_O_H(IIJB:IIJE,1:D%NKT)-1.))**0.58 + 0.37, 0.)
END WHERE
+!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
+!$mnh_expand_array(JIJ=IIJB:IIJE)
DO JK=IKTB+1,IKTE-1
- PMWTH(:,:,JK) = PMWTH(:,:,JK) * ZWSTAR(:,:)**2*ZTSTAR(:,:)
+ PMWTH(IIJB:IIJE,JK) = PMWTH(IIJB:IIJE,JK) * ZWSTAR(IIJB:IIJE)**2*ZTSTAR(IIJB:IIJE)
END DO
-PMWTH(:,:,IKE) = PMWTH(:,:,IKE) * ZWSTAR(:,:)**2*ZTSTAR(:,:)
-PMWTH(:,:,KKU) = PMWTH(:,:,KKU) * ZWSTAR(:,:)**2*ZTSTAR(:,:)
+PMWTH(IIJB:IIJE,IKE) = PMWTH(IIJB:IIJE,IKE) * ZWSTAR(IIJB:IIJE)**2*ZTSTAR(IIJB:IIJE)
+PMWTH(IIJB:IIJE,D%NKU) = PMWTH(IIJB:IIJE,D%NKU) * ZWSTAR(IIJB:IIJE)**2*ZTSTAR(IIJB:IIJE)
+!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
!----------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('TM06',1,ZHOOK_HANDLE)
diff --git a/src/common/turb/mode_update_lm.F90 b/src/common/turb/mode_update_lm.F90
index bbe5c9403e95ef29a89b99ae70cd2c7d2d424e59..5f3509398aa633335d26156660f8107ccbf351e1 100644
--- a/src/common/turb/mode_update_lm.F90
+++ b/src/common/turb/mode_update_lm.F90
@@ -6,7 +6,7 @@
MODULE MODE_UPDATE_LM
IMPLICIT NONE
CONTAINS
-SUBROUTINE UPDATE_LM(HLBCX,HLBCY,PLM,PLEPS)
+SUBROUTINE UPDATE_LM(D,HLBCX,HLBCY,PLM,PLEPS)
! #################################################################
!
!!**** *UPDATE_LM* - routine to set external points for mixing length
@@ -42,6 +42,7 @@ SUBROUTINE UPDATE_LM(HLBCX,HLBCY,PLM,PLEPS)
!* 0. DECLARATIONS
!
USE MODD_PARAMETERS
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
!
USE MODE_ll
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
@@ -51,11 +52,12 @@ IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X boundary type
CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y boundary type
!
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM ! mixing length
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS ! dissipative length
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PLM ! mixing length
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PLEPS ! dissipative length
!
!* 0.2 declarations of local variables
!
@@ -72,7 +74,10 @@ INTEGER :: IINFO_ll ! return code of parallel routine
!
!* 1. COMPUTE DIMENSIONS OF ARRAYS :
! ----------------------------
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IIB = D%NIB
+IIE = D%NIE
+IJB = D%NJB
+IJE = D%NJE
NULLIFY(TZLM_ll)
!
!-------------------------------------------------------------------------------
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index cb310eadfd0e35bae84296f84e2fef8e18c01f05..816b39604178c4e816b669f69ed9e82246ee86d0 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -246,7 +246,6 @@ USE MODD_BUDGET, ONLY: NBUDGET_U, NBUDGET_V, NBUDGET_W, NBUDGET_TH, NBUD
TBUDGETDATA, TBUDGETCONF_t
USE MODD_FIELD, ONLY: TFIELDDATA,TYPEREAL
USE MODD_IO, ONLY: TFILEDATA
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll
!
USE MODD_LES
USE MODD_IBM_PARAM_n, ONLY: LIBM, XIBM_LS, XIBM_XMUT
@@ -255,7 +254,7 @@ USE MODD_TURB_n, ONLY: TURB_t
!
USE MODE_BL89, ONLY: BL89
USE MODE_TURB_VER, ONLY : TURB_VER
-USE MODE_ROTATE_WIND, ONLY: ROTATE_WIND
+USE MODE_ROTATE_WIND, ONLY: ROTATE_WIND, UPDATE_ROTATE_WIND
USE MODE_TURB_HOR_SPLT, ONLY: TURB_HOR_SPLT
USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
USE MODE_RMC01, ONLY: RMC01
@@ -263,18 +262,16 @@ USE MODE_TM06, ONLY: TM06
USE MODE_UPDATE_LM, ONLY: UPDATE_LM
USE MODE_BUDGET, ONLY: BUDGET_STORE_INIT, BUDGET_STORE_END
USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
-USE MODE_ll, ONLY: ADD2DFIELD_ll, UPDATE_HALO_ll, CLEANLIST_ll, &
- LWEST_ll, LEAST_ll, LSOUTH_ll, LNORTH_ll
-USE MODE_SBL, ONLY: LMO
-USE MODE_SOURCES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT
+USE MODE_SBL_PHY, ONLY: LMO
+USE MODE_SOURCES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT_PHY
USE MODE_EMOIST, ONLY: EMOIST
USE MODE_ETHETA, ONLY: ETHETA
+USE MODE_IBM_MIXINGLENGTH, ONLY: IBM_MIXINGLENGTH
!
USE MODI_GRADIENT_W
USE MODI_GRADIENT_M
USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
-USE MODI_IBM_MIXINGLENGTH
USE MODI_LES_MEAN_SUBGRID
USE MODI_SHUMAN, ONLY : MZF, MXF, MYF
USE SHUMAN_PHY, ONLY : MZF_PHY
@@ -286,7 +283,7 @@ IMPLICIT NONE
!
!
!
-TYPE(DIMPHYEX_t), INTENT(IN) :: D
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
TYPE(CST_t), INTENT(IN) :: CST
TYPE(CSTURB_t), INTENT(IN) :: CSTURB
TYPE(TBUDGETCONF_t), INTENT(IN) :: BUCONF
@@ -503,7 +500,6 @@ REAL :: ZALPHA ! work coefficient :
!
REAL :: ZTIME1, ZTIME2
TYPE(TFIELDDATA) :: TZFIELD
-TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange (for UPDATE_ROTATE_WIND)
!
!* 1.PRELIMINARIES
! -------------
@@ -808,11 +804,11 @@ IF (ORMC01) THEN
ZUSTAR(:,:)=(PSFU(:,:)**2+PSFV(:,:)**2)**(0.25)
!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
IF (KRR>0) THEN
- CALL LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV,ZLMO)
+ CALL LMO(D,CST,ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV,ZLMO)
ELSE
ZRVM(:,:)=0.
ZSFRV(:,:)=0.
- CALL LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRVM,PSFTH,ZSFRV,ZLMO)
+ CALL LMO(D,CST,ZUSTAR,ZTHLM(:,:,IKB),ZRVM,PSFTH,ZSFRV,ZLMO)
END IF
CALL RMC01(D,CST,CSTURB,HTURBLEN,PZZ,PDXX,PDYY,PDZZ,PDIRCOSZW,PSBL_DEPTH,ZLMO,ZLM,ZLEPS)
END IF
@@ -828,14 +824,14 @@ END IF
! ----------------------------------------------------------
!
IF (HTURBDIM=="3DIM") THEN
- CALL UPDATE_LM(HLBCX,HLBCY,ZLM,ZLEPS)
+ CALL UPDATE_LM(D,HLBCX,HLBCY,ZLM,ZLEPS)
END IF
!
!* 3.9 Mixing length correction if immersed walls
! ------------------------------------------
!
IF (LIBM) THEN
- CALL IBM_MIXINGLENGTH(ZLM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET)
+ CALL IBM_MIXINGLENGTH(D,ZLM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET)
ENDIF
!----------------------------------------------------------------------------
!
@@ -848,13 +844,13 @@ ENDIF
!
!
IF (HPROGRAM/='AROME ') THEN
- CALL ROTATE_WIND(PUT,PVT,PWT, &
+ CALL ROTATE_WIND(D,PUT,PVT,PWT, &
PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
PCOSSLOPE,PSINSLOPE, &
PDXX,PDYY,PDZZ, &
ZUSLOPE,ZVSLOPE )
!
- CALL UPDATE_ROTATE_WIND(ZUSLOPE,ZVSLOPE)
+ CALL UPDATE_ROTATE_WIND(D,ZUSLOPE,ZVSLOPE,HLBCX,HLBCY)
ELSE
ZUSLOPE=PUT(IIB:IIE,IJB:IJE,D%NKA)
ZVSLOPE=PVT(IIB:IIE,IJB:IJE,D%NKA)
@@ -894,7 +890,7 @@ ZMR2(:,:,:) = 0. ! w'r'2
ZMTHR(:,:,:) = 0. ! w'th'r'
!
IF (HTOM=='TM06') THEN
- CALL TM06(D%NKA,D%NKU,D%NKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
+ CALL TM06(D,CST,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
!
ZFWTH = -GZ_M_W(D%NKA,D%NKU,D%NKL,ZMWTH,PDZZ) ! -d(w'2th' )/dz
!ZFWR = -GZ_M_W(D%NKA,D%NKU,D%NKL,ZMWR, PDZZ) ! -d(w'2r' )/dz
@@ -1273,7 +1269,7 @@ IF ( KRRL >= 1 ) THEN
END IF
! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
-CALL SOURCES_NEG_CORRECT(HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
+CALL SOURCES_NEG_CORRECT_PHY(D,KSV,HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
!----------------------------------------------------------------------------
!
!* 9. LES averaged surface fluxes
@@ -1345,73 +1341,6 @@ IF(PRESENT(PLEM)) PLEM(IIB:IIE,IJB:IJE,IKTB:IKTE) = ZLM(IIB:IIE,IJB:IJE,IKTB:IKT
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
CONTAINS
!
-!
-! ##############################################
- SUBROUTINE UPDATE_ROTATE_WIND(PUSLOPE,PVSLOPE)
-! ##############################################
-!!
-!!**** *UPDATE_ROTATE_WIND* routine to set rotate wind values at the border
-!
-!! AUTHOR
-!! ------
-!!
-!! P Jabouille *CNRM METEO-FRANCE
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 24/06/99
-!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
-!!
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PUSLOPE,PVSLOPE
-! tangential surface fluxes in the axes following the orography
-!
-IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',0,ZHOOK_HANDLE2)
-!
-!* 1 PROLOGUE
-!
-NULLIFY(TZFIELDS_ll)
-!
-! 2 Update halo if necessary
-!
-!!$IF (NHALO == 1) THEN
- CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' )
- CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' )
- CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_ll)
-!!$ENDIF
-!
-! 3 Boundary conditions for non cyclic case
-!
-IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
- PUSLOPE(D%NIB-1,:)=PUSLOPE(D%NIB,:)
- PVSLOPE(D%NIB-1,:)=PVSLOPE(D%NIB,:)
-END IF
-IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
- PUSLOPE(D%NIE+1,:)=PUSLOPE(D%NIE,:)
- PVSLOPE(D%NIE+1,:)=PVSLOPE(D%NIE,:)
-END IF
-IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
- PUSLOPE(:,D%NJB-1)=PUSLOPE(:,D%NJB)
- PVSLOPE(:,D%NJB-1)=PVSLOPE(:,D%NJB)
-END IF
-IF( HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
- PUSLOPE(:,D%NJE+1)=PUSLOPE(:,D%NJE)
- PVSLOPE(:,D%NJE+1)=PVSLOPE(:,D%NJE)
-END IF
-!
-IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',1,ZHOOK_HANDLE2)
-!
-END SUBROUTINE UPDATE_ROTATE_WIND
-!
! ########################################################################
SUBROUTINE COMPUTE_FUNCTION_THERMO(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
PLOCPEXN,PAMOIST,PATHETA )
diff --git a/src/mesonh/aux/sources_neg_correct.f90 b/src/mesonh/aux/sources_neg_correct.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9d2cffd55e23449660e63739b6b4f7d498dbfb45
--- /dev/null
+++ b/src/mesonh/aux/sources_neg_correct.f90
@@ -0,0 +1,383 @@
+!MNH_LIC Copyright 2020-2021 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.
+!-----------------------------------------------------------------
+! Author: P. Wautelet 25/06/2020 (deduplication of code from advection_metsv, resolved_cloud and turb)
+! Modifications:
+! P. Wautelet 30/06/2020: remove non-local corrections in resolved_cloud for NEGA => new local corrections here
+! J. Escobar 21/07/2020: bug <-> array of size(:,:,:,0) => return if krr=0
+! P. Wautelet 10/02/2021: budgets: add missing sources for NSV_C2R2BEG+3 budget
+!-----------------------------------------------------------------
+module mode_sources_neg_correct
+
+implicit none
+
+private
+
+public :: Sources_neg_correct,Sources_neg_correct_phy
+
+contains
+
+subroutine Sources_neg_correct_phy(D, KSV, hcloud, hbudname, KRR, ptstep, ppabst, ptht, prt, prths, prrs, prsvs, prhodj)
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KSV
+character(len=*), intent(in) :: hcloud ! Kind of cloud parameterization
+character(len=*), intent(in) :: hbudname ! Budget name
+integer, intent(in) :: KRR ! Number of moist variables
+real, intent(in) :: ptstep ! Timestep
+real, dimension(D%NIT,D%NJT,D%NKT), intent(in) :: ppabst ! Absolute pressure at time t
+real, dimension(D%NIT,D%NJT,D%NKT), intent(in) :: ptht ! Theta at time t
+real, dimension(D%NIT,D%NJT,D%NKT, KRR), intent(in) :: prt ! Moist variables at time t
+real, dimension(D%NIT,D%NJT,D%NKT), intent(inout) :: prths ! Source terms
+real, dimension(D%NIT,D%NJT,D%NKT, KRR), intent(inout) :: prrs ! Source terms
+real, dimension(D%NIT,D%NJT,D%NKT, KSV), intent(inout) :: prsvs ! Source terms
+real, dimension(D%NIT,D%NJT,D%NKT), intent(in), optional :: prhodj ! Dry density * jacobian
+!
+CALL SOURCES_NEG_CORRECT(HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHT,PRT,PRTHS,PRRS,PRSVS)
+!
+end subroutine Sources_neg_correct_phy
+!
+subroutine Sources_neg_correct( hcloud, hbudname, krr, ptstep, ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
+
+use modd_budget, only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_rr, lbudget_ri, &
+ lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv, &
+ NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RR, NBUDGET_RI, &
+ NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
+ tbudgets
+use modd_cst, only: xci, xcl, xcpd, xcpv, xlstt, xlvtt, xp00, xrd, xtt
+use modd_nsv, only: nsv_c2r2beg, nsv_c2r2end, nsv_lima_beg, nsv_lima_end, nsv_lima_nc, nsv_lima_nr, nsv_lima_ni
+use modd_param_lima, only: lcold_lima => lcold, lrain_lima => lrain, lspro_lima => lspro, lwarm_lima => lwarm, &
+ xctmin_lima => xctmin, xrtmin_lima => xrtmin
+
+use mode_budget, only: Budget_store_init, Budget_store_end
+use mode_msg
+
+implicit none
+
+character(len=*), intent(in) :: hcloud ! Kind of cloud parameterization
+character(len=*), intent(in) :: hbudname ! Budget name
+integer, intent(in) :: krr ! Number of moist variables
+real, intent(in) :: ptstep ! Timestep
+real, dimension(:, :, :), intent(in) :: ppabst ! Absolute pressure at time t
+real, dimension(:, :, :), intent(in) :: ptht ! Theta at time t
+real, dimension(:, :, :, :), intent(in) :: prt ! Moist variables at time t
+real, dimension(:, :, :), intent(inout) :: prths ! Source terms
+real, dimension(:, :, :, :), intent(inout) :: prrs ! Source terms
+real, dimension(:, :, :, :), intent(inout) :: prsvs ! Source terms
+real, dimension(:, :, :), intent(in), optional :: prhodj ! Dry density * jacobian
+
+integer :: ji, jj, jk
+integer :: jr
+integer :: jrmax
+integer :: jsv
+integer :: isv_lima_end
+real, dimension(:, :, :), allocatable :: zt, zexn, zlv, zls, zcph, zcor
+
+if ( krr == 0 ) return
+
+if ( hbudname /= 'NEADV' .and. hbudname /= 'NECON' .and. hbudname /= 'NEGA' .and. hbudname /= 'NETUR' ) &
+ call Print_msg( NVERB_WARNING, 'GEN', 'Sources_neg_correct', 'budget '//hbudname//' not yet tested' )
+
+if ( hcloud == 'LIMA' ) then
+ ! The negativity correction does not apply to the SPRO (supersaturation) variable which may be naturally negative
+ if ( lspro_lima ) then
+ isv_lima_end = nsv_lima_end - 1
+ else
+ isv_lima_end = nsv_lima_end
+ end if
+end if
+
+if ( hbudname /= 'NECON' .and. hbudname /= 'NEGA' ) then
+ if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+ hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+ if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
+ if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
+ if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
+ if ( lbudget_rr .and. &
+ ( hbudname /= 'NETUR' .or. &
+ ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
+ call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
+ IF (lbudget_ri .and. &
+ ( hbudname /= 'NETUR' .or. &
+ ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
+ call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
+ if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
+ if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
+ if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
+ end if
+
+ if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+ do ji = nsv_c2r2beg, nsv_c2r2end
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+ end do
+ end if
+ if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+ do ji = nsv_lima_beg, isv_lima_end
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+ end do
+ end if
+else !NECON + NEGA
+ if ( .not. present( prhodj ) ) &
+ call Print_msg( NVERB_FATAL, 'GEN', 'Sources_neg_correct', 'optional argument prhodj not present' )
+
+ if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+ hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+ if ( lbudget_th) call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) * prhodj(:, :, :) )
+ if ( lbudget_rv) call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) * prhodj(:, :, :) )
+ if ( lbudget_rc) call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) * prhodj(:, :, :) )
+ if ( lbudget_rr) call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) * prhodj(:, :, :) )
+ if ( lbudget_ri) call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) * prhodj(:, :, :) )
+ if ( lbudget_rs) call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) * prhodj(:, :, :) )
+ if ( lbudget_rg) call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) * prhodj(:, :, :) )
+ if ( lbudget_rh) call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) * prhodj(:, :, :) )
+ end if
+
+ if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+ do ji = nsv_c2r2beg, nsv_c2r2end
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+ end do
+ end if
+ if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+ do ji = nsv_lima_beg, isv_lima_end
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+ end do
+ end if
+end if
+
+allocate( zt ( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+allocate( zexn( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+allocate( zlv ( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+allocate( zcph( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+
+zexn(:, :, :) = ( ppabst(:, :, :) / xp00 ) ** (xrd / xcpd )
+zt (:, :, :) = ptht(:, :, :) * zexn(:, :, :)
+zlv (:, :, :) = xlvtt + ( xcpv - xcl ) * ( zt(:, :, :) - xtt )
+if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
+ allocate( zls( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+ zls(:, :, :) = xlstt + ( xcpv - xci ) * ( zt(:, :, :) - xtt )
+end if
+zcph(:, :, :) = xcpd + xcpv * prt(:, :, :, 1)
+
+deallocate( zt )
+
+CLOUD: select case ( hcloud )
+ case ( 'KESS' )
+ jrmax = Size( prrs, 4 )
+ do jr = 2, jrmax
+ where ( prrs(:, :, :, jr) < 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, jr) = 0.
+ end where
+ end do
+
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ end where
+
+
+ case( 'ICE3', 'ICE4' )
+ if ( hbudname == 'NETUR' ) then
+ jrmax = 4
+ else
+ jrmax = Size( prrs, 4 )
+ end if
+ do jr = 4, jrmax
+ where ( prrs(:, :, :, jr) < 0.)
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, jr) = 0.
+ end where
+ end do
+!
+! cloud
+ if ( hbudname == 'NETUR' ) then
+ jrmax = 2
+ else
+ jrmax = 3
+ end if
+ do jr = 2, jrmax
+ where ( prrs(:, :, :, jr) < 0.)
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, jr) = 0.
+ end where
+ end do
+!
+! if rc or ri are positive, we can correct negative rv
+! cloud
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ end where
+! ice
+ if ( krr > 3 ) then
+ allocate( zcor( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
+ zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
+ prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
+ end where
+ end if
+!
+!
+ case( 'C2R2', 'KHKO' )
+ where ( prrs(:, :, :, 2) < 0. .or. prsvs(:, :, :, nsv_c2r2beg + 1) < 0. )
+ prsvs(:, :, :, nsv_c2r2beg) = 0.
+ end where
+ do jsv = 2, 3
+ where ( prrs(:, :, :, jsv) < 0. .or. prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) < 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jsv)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jsv) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, jsv) = 0.
+ prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) = 0.
+ end where
+ end do
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ prsvs(:, :, :, nsv_c2r2beg + 1) = 0.
+ end where
+!
+!
+ case( 'LIMA' )
+! Correction where rc<0 or Nc<0
+ if ( lwarm_lima ) then
+ where ( prrs(:, :, :, 2) < xrtmin_lima(2) / ptstep .or. prsvs(:, :, :, nsv_lima_nc) < xctmin_lima(2) / ptstep )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ prsvs(:, :, :, nsv_lima_nc) = 0.
+ end where
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ prsvs(:, :, :, nsv_lima_nc) = 0.
+ end where
+ end if
+! Correction where rr<0 or Nr<0
+ if ( lwarm_lima .and. lrain_lima ) then
+ where ( prrs(:, :, :, 3) < xrtmin_lima(3) / ptstep .or. prsvs(:, :, :, nsv_lima_nr) < xctmin_lima(3) / ptstep )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 3)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 3) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 3) = 0.
+ prsvs(:, :, :, nsv_lima_nr) = 0.
+ end where
+ end if
+! Correction where ri<0 or Ni<0
+ if ( lcold_lima ) then
+ where ( prrs(:, :, :, 4) < xrtmin_lima(4) / ptstep .or. prsvs(:, :, :, nsv_lima_ni) < xctmin_lima(4) / ptstep )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 4)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 4) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 4) = 0.
+ prsvs(:, :, :, nsv_lima_ni) = 0.
+ end where
+ if ( hbudname /= 'NETUR' ) then
+ do jr = 5, Size( prrs, 4 )
+ where ( prrs(:, :, :, jr) < 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, jr) = 0.
+ end where
+ end do
+ end if
+ if(krr > 3) then
+ allocate( zcor( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
+ zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
+ prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
+ end where
+ deallocate( zcor )
+ end if
+ end if
+
+ prsvs(:, :, :, nsv_lima_beg : isv_lima_end) = Max( 0.0, prsvs(:, :, :, nsv_lima_beg : isv_lima_end) )
+
+end select CLOUD
+
+
+if ( hbudname /= 'NECON' .and. hbudname /= 'NEGA' ) then
+ if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+ hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+ if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
+ if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
+ if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
+ if ( lbudget_rr .and. &
+ ( hbudname /= 'NETUR' .or. &
+ ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
+ call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
+ IF (lbudget_ri .and. &
+ ( hbudname /= 'NETUR' .or. &
+ ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
+ call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
+ if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
+ if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
+ if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
+ end if
+
+ if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+ do ji = nsv_c2r2beg, nsv_c2r2end
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+ end do
+ end if
+ if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+ do ji = nsv_lima_beg, isv_lima_end
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+ end do
+ end if
+else !NECON + NEGA
+ if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+ hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+ if ( lbudget_th) call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) * prhodj(:, :, :) )
+ if ( lbudget_rv) call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) * prhodj(:, :, :) )
+ if ( lbudget_rc) call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) * prhodj(:, :, :) )
+ if ( lbudget_rr) call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) * prhodj(:, :, :) )
+ if ( lbudget_ri) call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) * prhodj(:, :, :) )
+ if ( lbudget_rs) call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) * prhodj(:, :, :) )
+ if ( lbudget_rg) call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) * prhodj(:, :, :) )
+ if ( lbudget_rh) call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) * prhodj(:, :, :) )
+ end if
+
+ if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+ do ji = nsv_c2r2beg, nsv_c2r2end
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+ end do
+ end if
+ if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+ do ji = nsv_lima_beg, isv_lima_end
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+ end do
+ end if
+end if
+
+end subroutine Sources_neg_correct
+
+end module mode_sources_neg_correct
diff --git a/src/mesonh/ext/ground_paramn.f90 b/src/mesonh/ext/ground_paramn.f90
index c213d7c6e8c8c981f3a005e5ae3bdbfc4552899e..1afd6b4d47b2bd2b2faeb9e118969f59ade6d9da 100644
--- a/src/mesonh/ext/ground_paramn.f90
+++ b/src/mesonh/ext/ground_paramn.f90
@@ -9,12 +9,14 @@ MODULE MODI_GROUND_PARAM_n
!
INTERFACE
!
- SUBROUTINE GROUND_PARAM_n( PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
+ SUBROUTINE GROUND_PARAM_n(D, PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
!* surface fluxes
! --------------
!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! surface flux of potential temperature (Km/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! surface flux of water vapor (m/s*kg/kg)
REAL, DIMENSION(:,:,:),INTENT(OUT):: PSFSV ! surface flux of scalar (m/s*kg/kg)
@@ -38,7 +40,7 @@ END INTERFACE
END MODULE MODI_GROUND_PARAM_n
!
! ######################################################################
- SUBROUTINE GROUND_PARAM_n( PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
+ SUBROUTINE GROUND_PARAM_n(D, PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
! #######################################################################
!
@@ -127,6 +129,7 @@ USE MODD_DYN, ONLY : XSEGLEN
!
USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_CST, ONLY : XP00, XCPD, XRD, XRV,XRHOLW, XDAY, XPI, XLVTT, XMD, XAVOGADRO
+USE MODD_DIMPHYEX, ONLY : DIMPHYEX_t
USE MODD_PARAMETERS, ONLY : JPVEXT, XUNDEF
USE MODD_DYN_n, ONLY : XTSTEP
USE MODD_CH_MNHC_n, ONLY : LUSECHEM
@@ -192,6 +195,7 @@ IMPLICIT NONE
!* surface fluxes
! --------------
!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! surface flux of potential temperature (Km/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! surface flux of water vapor (m/s*kg/kg)
REAL, DIMENSION(:,:,:),INTENT(OUT):: PSFSV ! surface flux of scalar (m/s*kg/kg)
@@ -417,7 +421,7 @@ END IF
! 1.3 Rotate the wind
! ---------------
!
-CALL ROTATE_WIND(XUT,XVT,XWT, &
+CALL ROTATE_WIND(D,XUT,XVT,XWT, &
XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
XCOSSLOPE,XSINSLOPE, &
XDXX,XDYY,XDZZ, &
diff --git a/src/mesonh/ext/phys_paramn.f90 b/src/mesonh/ext/phys_paramn.f90
index fe030a5192d348821cb0910d2708a9ed208d1523..e8ceabdd6e281809515126f69e4fd0a60f015cc1 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -1242,7 +1242,7 @@ IF (CSURF=='EXTE') THEN
DEALLOCATE( ZSAVE_INPRC,ZSAVE_PRCONV,ZSAVE_PRSCONV)
DEALLOCATE( ZSAVE_DIRFLASWD,ZSAVE_SCAFLASWD,ZSAVE_DIRSRFSWD)
END IF
- CALL GROUND_PARAM_n(ZSFTH, ZSFRV, ZSFSV, ZSFCO2, ZSFU, ZSFV, &
+ CALL GROUND_PARAM_n(YLDIMPHYEX,ZSFTH, ZSFRV, ZSFSV, ZSFCO2, ZSFU, ZSFV, &
ZDIR_ALB, ZSCA_ALB, ZEMIS, ZTSRAD )
!
IF (LIBM) THEN