From 89140b6d8a07e6c12908b58bf0beaf6f3f560ced Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Sat, 26 Feb 2022 00:46:35 +0100
Subject: [PATCH] Quentin 26/02/2022: missing bugfixes from last commit : 3 bug
bit-repro (add ifdef REPRO48/55) to be validated in compute_entr_detr and
compute_updraft
---
src/common/turb/mode_compute_entr_detr.F90 | 4 +
src/common/turb/mode_compute_updraft.F90 | 10 +-
src/mesonh/turb/mode_tridiag_massflux.f90 | 280 ---------------------
3 files changed, 13 insertions(+), 281 deletions(-)
delete mode 100644 src/mesonh/turb/mode_tridiag_massflux.f90
diff --git a/src/common/turb/mode_compute_entr_detr.F90 b/src/common/turb/mode_compute_entr_detr.F90
index 35bf18b95..7292dbb27 100644
--- a/src/common/turb/mode_compute_entr_detr.F90
+++ b/src/common/turb/mode_compute_entr_detr.F90
@@ -349,7 +349,11 @@ DO JLOOP=1,SIZE(OTEST)
(PRTM(JLOOP,KK)-ZDZ*(PRTM(JLOOP,KK)-PRTM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
(1. - ZKIC_INIT)*PRT_UP(JLOOP)
ELSE
+#ifdef REPRO55
+ ZMIXTHL(JLOOP) = 0.1
+#else
ZMIXTHL(JLOOP) = 300.
+#endif
ZMIXRT(JLOOP) = 0.1
ENDIF
ENDDO
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index cc5c9672a..3ad4b5584 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -319,7 +319,11 @@ IF (OENTR_DETR) THEN
ZSHEAR = 0. !no shear in bl89 mixing length
END IF
!
+#ifdef REPRO48
CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
+#else
+ CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.FALSE.,ZSHEAR,ZLUP)
+#endif
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
@@ -425,11 +429,15 @@ DO JK=KKB,KKE-KKL,KKL
ZMIX2(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*PENTR(JLOOP,JK) !&
ZMIX3_CLD(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZDETR_CLD(JLOOP,JK) !&
ZMIX2_CLD(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZENTR_CLD(JLOOP,JK)
-
+#ifdef REPRO48
PTHL_UP(JLOOP,JK+KKL)=(PTHL_UP(JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*ZMIX2(JLOOP)) &
/(1.+0.5*ZMIX2(JLOOP))
PRT_UP(JLOOP,JK+KKL) =(PRT_UP (JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PRTM(JLOOP,JK)*ZMIX2(JLOOP)) &
/(1.+0.5*ZMIX2(JLOOP))
+#else
+ PTHL_UP(JLOOP,JK+KKL)=PTHL_UP(JLOOP,JK)*EXP(-ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
+ PRT_UP(JLOOP,JK+KKL) =PRT_UP (JLOOP,JK)*EXP(-ZMIX2(JLOOP)) + PRTM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
+#endif
ENDIF
ENDDO
diff --git a/src/mesonh/turb/mode_tridiag_massflux.f90 b/src/mesonh/turb/mode_tridiag_massflux.f90
deleted file mode 100644
index 3e909c2c5..000000000
--- a/src/mesonh/turb/mode_tridiag_massflux.f90
+++ /dev/null
@@ -1,280 +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.
-MODULE MODE_TRIDIAG_MASSFLUX
-IMPLICIT NONE
-CONTAINS
-SUBROUTINE TRIDIAG_MASSFLUX(KKA,KKB,KKE,KKU,KKL,PVARM,PF,PDFDT,PTSTEP,PIMPL, &
- PDZZ,PRHODJ,PVARP )
-
- USE PARKIND1, ONLY : JPRB
- USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-! #################################################
-!
-!
-!!**** *TRIDIAG_MASSFLUX* - routine to solve a time implicit scheme
-!!
-!!
-!! PURPOSE
-!! -------
-! The purpose of this routine is to give a field PVARP at t+1, by
-! solving an implicit TRIDIAGonal system obtained by the
-! discretization of the vertical turbulent diffusion. It should be noted
-! that the degree of implicitness can be varied (PIMPL parameter) and that
-! the function of F(T) must have been linearized.
-! PVARP is localized at a mass point.
-!
-!!** METHOD
-!! ------
-!!
-!! [T(+) - T(-)]/2Dt = -d{ F + dF/dT *impl*[T(+) + T(-)] }/dz
-!!
-!! It is discretized as follows:
-!!
-!! PRHODJ(k)*PVARP(k)/PTSTEP
-!! =
-!! PRHODJ(k)*PVARM(k)/PTSTEP
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k+1)/PDZZ(k+1)
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARP(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k) /PDZZ(k+1)
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARP(k) /PDZZ(k+1)
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARM(k) /PDZZ(k)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARP(k) /PDZZ(k)
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARM(k-1)/PDZZ(k)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARP(k-1)/PDZZ(k)
-!!
-!!
-!! The system to solve is:
-!!
-!! A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k)
-!!
-!!
-!! The RHS of the linear system in PVARP writes:
-!!
-!! y(k) = PRHODJ(k)*PVARM(k)/PTSTEP
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k) /PDZZ(k+1)
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARM(k) /PDZZ(k)
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) * PVARM(k-1)/PDZZ(k)
-!!
-!!
-!! Then, the classical TRIDIAGonal algorithm is used to invert the
-!! implicit operator. Its matrix is given by:
-!!
-!! ( b(KKB) c(KKB) 0 0 0 0 0 0 )
-!! ( a(KKB+1) b(KKB+1) c(KKB+1) 0 ... 0 0 0 0 )
-!! ( 0 a(KKB+2) b(KKB+2) c(KKB+2). 0 0 0 0 )
-!! .......................................................................
-!! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 )
-!! .......................................................................
-!! ( 0 0 0 0 0 ...a(KKE-1) b(KKE-1) c(KKE-1))
-!! ( 0 0 0 0 0 ... 0 a(KKE) b(KKE) )
-!!
-!! KKB and KKE represent the first and the last inner mass levels of the
-!! model. The coefficients are:
-!!
-!! a(k) = - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) /PDZZ(k)
-!! b(k) = PRHODJ(k) / PTSTEP
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1)/PDZZ(k+1)
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k) /PDZZ(k)
-!! c(k) = + (PRHODJ(k+1)+PRHODJ(k) )/2. * 0.5*PIMPL* PDFDT(k+1)/PDZZ(k+1)
-!!
-!! for all k /= KKB or KKE
-!!
-!!
-!! b(KKB) = PRHODJ(KKB) / PTSTEP
-!! +(PRHODJ(KKB+1)+PRHODJ(KKB))/2.*0.5*PIMPL*PDFDT(KKB+1)/PDZZ(KKB+1)
-!! c(KKB) = +(PRHODJ(KKB+1)+PRHODJ(KKB))/2.*0.5*PIMPL*PDFDT(KKB+1)/PDZZ(KKB+1)
-!!
-!! b(KKE) = PRHODJ(KKE) / PTSTEP
-!! -(PRHODJ(KKE)+PRHODJ(KKE-1))/2.*0.5*PIMPL*PDFDT(KKE)/PDZZ(KKE)
-!! a(KKE) = -(PRHODJ(KKE)+PRHODJ(KKE-1))/2.*0.5*PIMPL*PDFDT(KKE)/PDZZ(KKE)
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! NONE
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!! Press et al: Numerical recipes (1986) Cambridge Univ. Press
-!!
-!! AUTHOR
-!! ------
-!! V. Masson and S. Malardel * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 07/2006
-!! V.Masson : Optimization
-!! S. Riette Jan 2012: support for both order of vertical levels
-!! ---------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-!
-USE MODD_PARAMETERS, ONLY: JPVEXT
-USE MODI_SHUMAN_MF
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-INTEGER, INTENT(IN) :: KKA ! near ground array index
-INTEGER, INTENT(IN) :: KKB ! near ground physical index
-INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
-INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
-INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point
-REAL, DIMENSION(:,:), INTENT(IN) :: PF ! flux in dT/dt=-dF/dz at flux point
-REAL, DIMENSION(:,:), INTENT(IN) :: PDFDT ! dF/dT at flux point
-REAL, INTENT(IN) :: PTSTEP ! Double time step
-REAL, INTENT(IN) :: PIMPL ! implicit weight
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! Dz at flux point
-REAL, DIMENSION(:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J at mass point
-!
-REAL, DIMENSION(:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point
-!
-!
-!* 0.2 declarations of local variables
-!
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZRHODJ_DFDT_O_DZ
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZMZM_RHODJ
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZA, ZB, ZC
-REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZY ,ZGAM
- ! RHS of the equation, 3D work array
-REAL, DIMENSION(SIZE(PVARM,1)) :: ZBET
- ! 2D work array
-INTEGER :: JK ! loop counter
-!
-! ---------------------------------------------------------------------------
-!
-!* 1. Preliminaries
-! -------------
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('TRIDIAG_MASSFLUX',0,ZHOOK_HANDLE)
-ZMZM_RHODJ = MZM_MF(KKA, KKU, KKL, PRHODJ)
-ZRHODJ_DFDT_O_DZ = ZMZM_RHODJ*PDFDT/PDZZ
-!
-ZA=0.
-ZB=0.
-ZC=0.
-ZY=0.
-!
-!
-!* 2. COMPUTE THE RIGHT HAND SIDE
-! ---------------------------
-!
-ZY(:,KKB) = PRHODJ(:,KKB)*PVARM(:,KKB)/PTSTEP &
- - ZMZM_RHODJ(:,KKB+KKL) * PF(:,KKB+KKL)/PDZZ(:,KKB+KKL) &
- + ZMZM_RHODJ(:,KKB ) * PF(:,KKB )/PDZZ(:,KKB ) &
- + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL * PVARM(:,KKB+KKL) &
- + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL * PVARM(:,KKB )
-!
-DO JK=2+JPVEXT,SIZE(ZY,2)-JPVEXT-1
- ZY(:,JK) = PRHODJ(:,JK)*PVARM(:,JK)/PTSTEP &
- - ZMZM_RHODJ(:,JK+KKL) * PF(:,JK+KKL)/PDZZ(:,JK+KKL) &
- + ZMZM_RHODJ(:,JK ) * PF(:,JK )/PDZZ(:,JK ) &
- + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL * PVARM(:,JK+KKL) &
- + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL * PVARM(:,JK ) &
- - ZRHODJ_DFDT_O_DZ(:,JK ) * 0.5*PIMPL * PVARM(:,JK ) &
- - ZRHODJ_DFDT_O_DZ(:,JK ) * 0.5*PIMPL * PVARM(:,JK-KKL)
-END DO
-!
-IF (JPVEXT==0) THEN
- ZY(:,KKE) = PRHODJ(:,KKE)*PVARM(:,KKE)/PTSTEP
-ELSE
- ZY(:,KKE) = PRHODJ(:,KKE)*PVARM(:,KKE)/PTSTEP &
- - ZMZM_RHODJ(:,KKE+KKL) * PF(:,KKE+KKL)/PDZZ(:,KKE+KKL) &
- + ZMZM_RHODJ(:,KKE ) * PF(:,KKE )/PDZZ(:,KKE ) &
- - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL * PVARM(:,KKE ) &
- - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL * PVARM(:,KKE-KKL)
-ENDIF
-!
-!
-!* 3. INVERSION OF THE TRIDIAGONAL SYSTEM
-! -----------------------------------
-!
-IF ( PIMPL > 1.E-10 ) THEN
-!
-!* 3.1 arrays A, B, C
-! --------------
-!
- ZB(:,KKB) = PRHODJ(:,KKB)/PTSTEP &
- + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL
- ZC(:,KKB) = ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL
-
- DO JK=2+JPVEXT,SIZE(ZY,2)-JPVEXT-1
- ZA(:,JK) = - ZRHODJ_DFDT_O_DZ(:,JK ) * 0.5*PIMPL
- ZB(:,JK) = PRHODJ(:,JK)/PTSTEP &
- + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL &
- - ZRHODJ_DFDT_O_DZ(:,JK ) * 0.5*PIMPL
- ZC(:,JK) = ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL
- END DO
-
- ZA(:,KKE) = - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL
- ZB(:,KKE) = PRHODJ(:,KKE)/PTSTEP &
- - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL
-!
-!* 3.2 going up
-! --------
-!
- ZBET(:) = ZB(:,KKB) ! bet = b(KKB)
- PVARP(:,KKB) = ZY(:,KKB) / ZBET(:)
-
- !
- DO JK = KKB+KKL,KKE-KKL,KKL
- ZGAM(:,JK) = ZC(:,JK-KKL) / ZBET(:)
- ! gam(k) = c(k-1) / bet
- ZBET(:) = ZB(:,JK) - ZA(:,JK) * ZGAM(:,JK)
- ! bet = b(k) - a(k)* gam(k)
- PVARP(:,JK)= ( ZY(:,JK) - ZA(:,JK) * PVARP(:,JK-KKL) ) / ZBET(:)
- ! res(k) = (y(k) -a(k)*res(k-1))/ bet
- END DO
- ! special treatment for the last level
- ZGAM(:,KKE) = ZC(:,KKE-KKL) / ZBET(:)
- ! gam(k) = c(k-1) / bet
- ZBET(:) = ZB(:,KKE) - ZA(:,KKE) * ZGAM(:,KKE)
- ! bet = b(k) - a(k)* gam(k)
- PVARP(:,KKE)= ( ZY(:,KKE) - ZA(:,KKE) * PVARP(:,KKE-KKL) ) / ZBET(:)
- ! res(k) = (y(k) -a(k)*res(k-1))/ bet
-!
-!* 3.3 going down
-! ----------
-!
- DO JK = KKE-KKL,KKB,-KKL
- PVARP(:,JK) = PVARP(:,JK) - ZGAM(:,JK+KKL) * PVARP(:,JK+KKL)
- END DO
-!
-!
-ELSE
- !!! EXPLICIT FORMULATION
- !
- DO JK=1+JPVEXT,SIZE(PVARP,2)-JPVEXT
- PVARP(:,JK) = ZY(:,JK) * PTSTEP / PRHODJ(:,JK)
- ENDDO
- !
-END IF
-!
-!
-!* 4. FILL THE UPPER AND LOWER EXTERNAL VALUES
-! ----------------------------------------
-!
-PVARP(:,KKA)=PVARP(:,KKB)
-PVARP(:,KKU)=PVARP(:,KKE)
-!
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('TRIDIAG_MASSFLUX',1,ZHOOK_HANDLE)
-END SUBROUTINE TRIDIAG_MASSFLUX
-END MODULE MODE_TRIDIAG_MASSFLUX
--
GitLab