From 9746760225e99c2ffb1fcbae2659b784bda0b3b2 Mon Sep 17 00:00:00 2001 From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr> Date: Tue, 30 Jun 2020 10:31:44 +0200 Subject: [PATCH] Philippe 30/06/2020: move removal of negative scalar variables to Sources_neg_correct --- src/MNH/advection_metsv.f90 | 186 +---------------- src/MNH/resolved_cloud.f90 | 350 +------------------------------- src/MNH/sources_neg_correct.f90 | 292 ++++++++++++++++++++++++++ src/MNH/turb.f90 | 184 +---------------- 4 files changed, 315 insertions(+), 697 deletions(-) create mode 100644 src/MNH/sources_neg_correct.f90 diff --git a/src/MNH/advection_metsv.f90 b/src/MNH/advection_metsv.f90 index 6362da54e..cf7cf6e9e 100644 --- a/src/MNH/advection_metsv.f90 +++ b/src/MNH/advection_metsv.f90 @@ -139,8 +139,8 @@ END MODULE MODI_ADVECTION_METSV !! 03/2020 (B.Vie) : LIMA negativity checks after turbulence, advection and !! microphysics budgets ! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices -! P. Wautelet + Benoît Vié 11/06/2020: improve removal of negative scalar variables + adapt the corresponding budgets -! P. Wautelet + Benoît Vié 23/06/2020: improve removal of negative scalar variables +! P. Wautelet + Benoît Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets +! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -166,6 +166,7 @@ USE MODE_FMWRIT USE MODE_IO_ll USE MODE_ll USE MODE_MSG +use mode_sources_neg_correct, only: Sources_neg_correct ! USE MODI_ADV_BOUNDARIES USE MODI_BUDGET @@ -692,185 +693,10 @@ IF (KRR>=7.AND.LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'ADV_BU_RRH') DO JSV=1,KSV IF (LBUDGET_SV) CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'ADV_BU_RSV') END DO -! -SELECT CASE ( HCLOUD ) - CASE('KESS') - ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - DO JR=2,SIZE(PRRS,4) - WHERE (PRRS(:,:,:,JR) < 0.) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,JR) - PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,JR) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRRS(:,:,:,JR) = 0.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') - ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - DO JR=4,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 - DO JR=2,3 - 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 - 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') - ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - WHERE (PRRS(:,:,:,2) < 0. .OR. PRSVS(:,:,:,NSV_C2R2BEG+1) < 0.) - PRSVS(:,:,:,NSV_C2R2BEG) = 0.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.0 - PRSVS(:,:,:,NSV_C2R2BEG-1+JSV) = 0.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') - ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) -! Correction where rc<0 or Nc<0 - IF (LWARM) THEN - WHERE (PRRS(:,:,:,2) < XRTMIN(2)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NC) < XCTMIN(2)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRRS(:,:,:,2) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NC) = 0.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.0 - END WHERE - END IF -! Correction where rr<0 or Nr<0 - IF (LWARM .AND. LRAIN) THEN - WHERE (PRRS(:,:,:,3) < XRTMIN(3)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NR) < XCTMIN(3)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,3) - PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,3) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRRS(:,:,:,3) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NR) = 0.0 - END WHERE - END IF -! Correction where ri<0 or Ni<0 - IF (LCOLD) THEN - WHERE (PRRS(:,:,:,4) < XRTMIN(4)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NI) < XCTMIN(4)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,4) - PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRRS(:,:,:,4) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NI) = 0.0 - END WHERE - 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 - IF(KRR > 3) THEN - 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 - END IF -! - PRSVS(:, :, :, NSV_LIMA_BEG:NSV_LIMA_END) = MAX( 0.0, PRSVS(:, :, :, NSV_LIMA_BEG:NSV_LIMA_END) ) -! -END SELECT -! -IF ( HCLOUD == 'KESS' .OR. HCLOUD == 'ICE3' .OR. HCLOUD == 'ICE4' .OR. & - HCLOUD == 'KHKO' .OR. HCLOUD == 'C2R2' .OR. HCLOUD == 'LIMA' ) THEN - IF (LBUDGET_TH) CALL BUDGET (PRTHS(:,:,:) , 4,'NEADV_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NEADV_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NEADV_BU_RRC') - IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3), 8,'NEADV_BU_RRR') - IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4), 9,'NEADV_BU_RRI') - IF (LBUDGET_RS) CALL BUDGET (PRRS(:,:,:,5),10,'NEADV_BU_RRS') - IF (LBUDGET_RG) CALL BUDGET (PRRS(:,:,:,6),11,'NEADV_BU_RRG') - IF (LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'NEADV_BU_RRH') -END IF -IF (LBUDGET_SV .AND. (HCLOUD == 'C2R2' .OR. HCLOUD == 'KHKO')) THEN - DO JI = 1, 3 - CALL BUDGET ( PRSVS(:,:,:,NSV_C2R2BEG-1+JI),12+NSV_C2R2BEG-1+JI,'NEADV_BU_RSV') - END DO -END IF -IF (LBUDGET_SV .AND. HCLOUD == 'LIMA') THEN - DO JI = NSV_LIMA_BEG, NSV_LIMA_END - CALL BUDGET ( PRSVS(:,:,:,JI),12+JI,'NEADV_BU_RSV') - END DO -END IF +! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets +call Sources_neg_correct( hcloud, 'NEADV', krr, ptstep, ppabst, ptht, prt, prths, prrs, prsvs ) + !------------------------------------------------------------------------------- ! END SUBROUTINE ADVECTION_METSV diff --git a/src/MNH/resolved_cloud.f90 b/src/MNH/resolved_cloud.f90 index 8ae441ee0..7b313bbb7 100644 --- a/src/MNH/resolved_cloud.f90 +++ b/src/MNH/resolved_cloud.f90 @@ -272,10 +272,9 @@ END MODULE MODI_RESOLVED_CLOUD !! B.Vié 03/03/2020 : use DTHRAD instead of dT/dt in Smax diagnostic computation ! P. Wautelet 11/06/2020: bugfix: correct ZSVS array indices ! P. Wautelet 11/06/2020: bugfix: add "Non local correction for precipitating species" for ICE4 -! P. Wautelet + Benoit Vié 11/06/2020: improve removal of negative scalar variables + adapt the corresponding budgets -! P. Wautelet + Benoit Vié 18/06/2020: improve removal of negative scalar variables: use ZEXN instead of PEXNREF +! P. Wautelet + Benoit Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets ! P. Wautelet 23/06/2020: remove ZSVS and ZSVT to improve code readability -! P. Wautelet + Benoît Vié 23/06/2020: improve removal of negative scalar variables +! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -299,10 +298,11 @@ USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN USE MODD_SALT, ONLY: LSALT ! USE MODE_ll +use mode_sources_neg_correct, only: Sources_neg_correct ! USE MODI_BUDGET USE MODI_C2R2_ADJUST -USE MODI_C3R5_ADJUST +! USE MODI_C3R5_ADJUST USE MODI_FAST_TERMS USE MODI_GET_HALO USE MODI_ICE_ADJUST @@ -571,17 +571,6 @@ IF (HCLOUD == 'C2R2' .OR. HCLOUD == 'C3R5' .OR. HCLOUD == 'KHKO' & PSVT(:,:,IKE+1,ISVBEG:ISVEND) = PSVT(:,:,IKE,ISVBEG:ISVEND) ENDIF ! -! personal comment: tranfering these variables to the -! microphysical routines would save -! computing time -! -ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) -ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:) -ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT) -ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT) -ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) -! -! !* 3. REMOVE NEGATIVE VALUES ! ---------------------- ! @@ -622,164 +611,8 @@ END IF ! !* 3.2 Adjustement for liquid and solid cloud ! -SELECT CASE ( HCLOUD ) - CASE('KESS') - WHERE (PRS(:,:,:,2) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0.0 - END WHERE - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - END WHERE -! -! -! - CASE('C2R2','KHKO') - WHERE (PRS(:,:,:,2) < 0. .OR. PSVS(:,:,:,NSV_C2R2BEG+1) < 0.) - PSVS(:,:,:,NSV_C2R2BEG+1) = 0.0 - END WHERE - DO JSV = 2, 3 - WHERE (PRS(:,:,:,JSV) < 0. .OR. PSVS(:,:,:,NSV_C2R2BEG-1+JSV) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JSV) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JSV) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JSV) = 0.0 - PSVS(:,:,:,NSV_C2R2BEG-1+JSV) = 0.0 - END WHERE - ENDDO - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - PSVS(:,:,:,NSV_C2R2BEG+1) = 0. - END WHERE -! -! -! - CASE('ICE3','ICE4') - WHERE (PRS(:,:,:,4) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,4) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = 0. - END WHERE -! -! cloud - WHERE (PRS(:,:,:,2) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - END WHERE -! -! if rc or ri are positive, we can correct negative rv -! cloud - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - END WHERE -! ice - IF(KRR > 3) THEN - WHERE ((PRS(:,:,:,1) < 0.).AND.(PRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRS(:,:,:,1),PRS(:,:,:,4)) - PRS(:,:,:,1) = PRS(:,:,:,1) + ZCOR(:,:,:) - PTHS(:,:,:) = PTHS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = PRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF -! - CASE('LIMA') -! Correction where rc<0 or Nc<0 - IF (OWARM) THEN - WHERE (PRS(:,:,:,2) < YRTMIN(2)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NC) < YCTMIN(2)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0.0 - PSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - PSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - END IF -! Correction where rr<0 or Nr<0 - IF (OWARM .AND. ORAIN) THEN - WHERE (PRS(:,:,:,3) < YRTMIN(3)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NR) < YCTMIN(3)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,3) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,3) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,3) = 0.0 - PSVS(:,:,:,NSV_LIMA_NR) = 0.0 - END WHERE - END IF -! Correction where ri<0 or Ni<0 - IF (LCOLD) THEN - WHERE (PRS(:,:,:,4) < YRTMIN(4)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NI) < YCTMIN(4)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,4) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = 0.0 - PSVS(:,:,:,NSV_LIMA_NI) = 0.0 - END WHERE - IF(KRR > 3) THEN - WHERE ((PRS(:,:,:,1) < 0.).AND.(PRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRS(:,:,:,1),PRS(:,:,:,4)) - PRS(:,:,:,1) = PRS(:,:,:,1) + ZCOR(:,:,:) - PTHS(:,:,:) = PTHS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = PRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF - END IF -! if rc or ri are positive, we can correct negative rv -! cloud -! ice -! - PSVS(:,:,:,NSV_LIMA_BEG:NSV_LIMA_END) = MAX( 0.0,PSVS(:,:,:,NSV_LIMA_BEG:NSV_LIMA_END) ) -! -END SELECT -! -! -!* 3.3 STORE THE BUDGET TERMS -! ---------------------- -! -IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:) * PRHODJ(:,:,:), 4,'NEGA_BU_RTH') -IF (LBUDGET_RV) CALL BUDGET (PRS(:,:,:,1) * PRHODJ(:,:,:), 6,'NEGA_BU_RRV') -IF (LBUDGET_RC) CALL BUDGET (PRS(:,:,:,2) * PRHODJ(:,:,:), 7,'NEGA_BU_RRC') -IF (LBUDGET_RR) CALL BUDGET (PRS(:,:,:,3) * PRHODJ(:,:,:), 8,'NEGA_BU_RRR') -IF (LBUDGET_RI) CALL BUDGET (PRS(:,:,:,4) * PRHODJ(:,:,:) ,9,'NEGA_BU_RRI') -IF (LBUDGET_RS) CALL BUDGET (PRS(:,:,:,5) * PRHODJ(:,:,:),10,'NEGA_BU_RRS') -IF (LBUDGET_RG) CALL BUDGET (PRS(:,:,:,6) * PRHODJ(:,:,:),11,'NEGA_BU_RRG') -IF (LBUDGET_RH) CALL BUDGET (PRS(:,:,:,7) * PRHODJ(:,:,:),12,'NEGA_BU_RRH') - -IF (LBUDGET_SV .AND. (HCLOUD == 'C2R2' .OR. HCLOUD == 'KHKO')) THEN - DO JSV = 1, 3 - CALL BUDGET ( PSVS(:,:,:,NSV_C2R2BEG-1+JSV),12+NSV_C2R2BEG-1+JSV,'NEGA_BU_RSV') - END DO -END IF -IF (LBUDGET_SV .AND. HCLOUD == 'C3R5') THEN - DO JSV = 1, 4 - CALL BUDGET ( PSVS(:,:,:,NSV_C2R2BEG-1+JSV),12+NSV_C2R2BEG-1+JSV,'NEGA_BU_RSV') - END DO -END IF -IF (LBUDGET_SV .AND. HCLOUD == 'LIMA') THEN - DO JSV = NSV_LIMA_BEG, NSV_LIMA_END - CALL BUDGET ( PSVS(:,:,:,JSV),12+JSV,'NEGA_BU_RSV') - END DO -END IF +! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets +call Sources_neg_correct( hcloud, 'NEGA', krr, ptstep, ppabst, ptht, prt, pths, prs, psvs, prhodj ) ! !* 3.4 Limitations of Na and Nc to the CCN max number concentration ! @@ -1129,175 +962,10 @@ IF(HCLOUD=='ICE3' .OR. HCLOUD=='ICE4' ) THEN ENDWHERE ENDIF ENDIF -! -! -SELECT CASE ( HCLOUD ) - CASE('KESS') - DO JRR=2,SIZE(PRS,4) - WHERE (PRS(:,:,:,JRR) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JRR) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JRR) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JRR) = 0.0 - END WHERE - END DO - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - END WHERE -! -! - CASE('ICE3','ICE4') - DO JRR=4,SIZE(PRS,4) - WHERE (PRS(:,:,:,JRR) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JRR) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JRR) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JRR) = 0. - END WHERE - END DO -! -! cloud - DO JRR=2,3 - WHERE (PRS(:,:,:,JRR) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JRR) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JRR) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JRR) = 0. - END WHERE - END DO -! -! if rc or ri are positive, we can correct negative rv -! cloud - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - END WHERE -! ice - IF(KRR > 3) THEN - WHERE ((PRS(:,:,:,1) < 0.).AND.(PRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRS(:,:,:,1),PRS(:,:,:,4)) - PRS(:,:,:,1) = PRS(:,:,:,1) + ZCOR(:,:,:) - PTHS(:,:,:) = PTHS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = PRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF -! -! - CASE('C2R2','KHKO') - WHERE (PRS(:,:,:,2) < 0. .OR. PSVS(:,:,:,NSV_C2R2BEG+1) < 0.) - PSVS(:,:,:,NSV_C2R2BEG) = 0.0 - END WHERE - DO JSV = 2, 3 - WHERE (PRS(:,:,:,JSV) < 0. .OR. PSVS(:,:,:,NSV_C2R2BEG-1+JSV) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JSV) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JSV) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JSV) = 0.0 - PSVS(:,:,:,NSV_C2R2BEG-1+JSV) = 0.0 - END WHERE - ENDDO - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - PSVS(:,:,:,NSV_C2R2BEG+1) = 0. - END WHERE -! - CASE('LIMA') -! Correction where rc<0 or Nc<0 - IF (OWARM) THEN - WHERE (PRS(:,:,:,2) < YRTMIN(2)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NC) < YCTMIN(2)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0.0 - PSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - WHERE ((PRS(:,:,:,1) <0.) .AND. (PRS(:,:,:,2)> 0.) ) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,2) = 0. - PSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - END IF -! Correction where rr<0 or Nr<0 - IF (OWARM .AND. ORAIN) THEN - WHERE (PRS(:,:,:,3) < YRTMIN(3)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NR) < YCTMIN(3)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,3) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,3) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,3) = 0.0 - PSVS(:,:,:,NSV_LIMA_NR) = 0.0 - END WHERE - END IF -! Correction where ri<0 or Ni<0 - IF (LCOLD) THEN - WHERE (PRS(:,:,:,4) < YRTMIN(4)/PTSTEP .OR. PSVS(:,:,:,NSV_LIMA_NI) < YCTMIN(4)/PTSTEP) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,4) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = 0.0 - PSVS(:,:,:,NSV_LIMA_NI) = 0.0 - END WHERE - DO JRR=5,SIZE(PRS,4) - WHERE (PRS(:,:,:,JRR) < 0.) - PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,JRR) - PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,JRR) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,JRR) = 0. - END WHERE - END DO - IF(KRR > 3) THEN - WHERE ((PRS(:,:,:,1) < 0.).AND.(PRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRS(:,:,:,1),PRS(:,:,:,4)) - PRS(:,:,:,1) = PRS(:,:,:,1) + ZCOR(:,:,:) - PTHS(:,:,:) = PTHS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXN(:,:,:) - PRS(:,:,:,4) = PRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF - END IF -! - PSVS(:,:,:,NSV_LIMA_BEG:NSV_LIMA_END) = MAX( 0.0,PSVS(:,:,:,NSV_LIMA_BEG:NSV_LIMA_END) ) -! -END SELECT -! -! -!* 3.3 STORE THE BUDGET TERMS -! ---------------------- -! -IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:) * PRHODJ(:,:,:), 4,'NECON_BU_RTH') -IF (LBUDGET_RV) CALL BUDGET (PRS(:,:,:,1) * PRHODJ(:,:,:), 6,'NECON_BU_RRV') -IF (LBUDGET_RC) CALL BUDGET (PRS(:,:,:,2) * PRHODJ(:,:,:), 7,'NECON_BU_RRC') -IF (LBUDGET_RR) CALL BUDGET (PRS(:,:,:,3) * PRHODJ(:,:,:), 8,'NECON_BU_RRR') -IF (LBUDGET_RI) CALL BUDGET (PRS(:,:,:,4) * PRHODJ(:,:,:) ,9,'NECON_BU_RRI') -IF (LBUDGET_RS) CALL BUDGET (PRS(:,:,:,5) * PRHODJ(:,:,:),10,'NECON_BU_RRS') -IF (LBUDGET_RG) CALL BUDGET (PRS(:,:,:,6) * PRHODJ(:,:,:),11,'NECON_BU_RRG') -IF (LBUDGET_RH) CALL BUDGET (PRS(:,:,:,7) * PRHODJ(:,:,:),12,'NECON_BU_RRH') -IF (LBUDGET_SV .AND. (HCLOUD == 'C2R2' .OR. HCLOUD == 'KHKO')) THEN - DO JSV = 1, 3 - CALL BUDGET ( PSVS(:,:,:,NSV_C2R2BEG-1+JSV),12+NSV_C2R2BEG-1+JSV,'NECON_BU_RSV') - END DO -END IF -IF (LBUDGET_SV .AND. HCLOUD == 'C3R5') THEN - DO JSV = 1, 4 - CALL BUDGET ( PSVS(:,:,:,NSV_C2R2BEG-1+JSV),12+NSV_C2R2BEG-1+JSV,'NECON_BU_RSV') - END DO -END IF -IF (LBUDGET_SV .AND. HCLOUD == 'LIMA') THEN - DO JSV = NSV_LIMA_BEG, NSV_LIMA_END - CALL BUDGET ( PSVS(:,:,:,JSV),12+JSV,'NECON_BU_RSV') - END DO -END IF +! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets +call Sources_neg_correct( hcloud, 'NECON', krr, ptstep, ppabst, ptht, prt, pths, prs, psvs, prhodj ) + !------------------------------------------------------------------------------- ! ! diff --git a/src/MNH/sources_neg_correct.f90 b/src/MNH/sources_neg_correct.f90 new file mode 100644 index 000000000..234ea49fe --- /dev/null +++ b/src/MNH/sources_neg_correct.f90 @@ -0,0 +1,292 @@ +!MNH_LIC Copyright 2020-2020 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: +! +!----------------------------------------------------------------- +module mode_sources_neg_correct + +implicit none + +private + +public :: Sources_neg_correct + +contains + +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 +use modd_cst, only: xci, xcl, xcpd, xcpv, xlstt, xlvtt, xp00, xrd, xtt +use modd_nsv, only: nsv_c2r2beg, 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, lwarm_lima => lwarm, & + xctmin_lima => xctmin, xrtmin_lima => xrtmin + +use mode_msg + +use modi_budget + +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 +real, dimension(:, :, :), allocatable :: zt, zexn, zlv, zls, zcph, zcor + +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' ) + +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' ) + if ( hbudname == 'NEGA' ) then + jrmax = 2 + else + jrmax = Size( prrs, 4 ) + 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 + + 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' .or. hbudname == 'NEGA' ) 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' .or. hbudname == 'NEGA' ) 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' .and. hbudname /= 'NEGA' ) 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 : nsv_lima_end) = Max( 0.0, prsvs(:, :, :, nsv_lima_beg : nsv_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( prths(:, :, :), 4, Trim( hbudname )//'_BU_RTH' ) + if ( lbudget_rv ) call Budget( prrs (:, :, :, 1), 6, Trim( hbudname )//'_BU_RRV' ) + if ( lbudget_rc ) call Budget( prrs (:, :, :, 2), 7, Trim( hbudname )//'_BU_RRC' ) + if ( lbudget_rr .and. & + ( hbudname /= 'NETUR' .or. & + ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) & + call Budget( prrs(:, :, :, 3), 8, Trim( hbudname )//'_BU_RRR' ) + IF (lbudget_ri .and. & + ( hbudname /= 'NETUR' .or. & + ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) & + call Budget( prrs(:, :, :, 4), 9, Trim( hbudname )//'_BU_RRI' ) + IF (lbudget_rs .and. hbudname /= 'NETUR' ) call Budget( prrs(:, :, :, 5), 10, Trim( hbudname )//'_BU_RRS' ) + IF (lbudget_rg .and. hbudname /= 'NETUR' ) call Budget( prrs(:, :, :, 6), 11, Trim( hbudname )//'_BU_RRG' ) + IF (lbudget_rh .and. hbudname /= 'NETUR' ) call Budget( prrs(:, :, :, 7), 12, Trim( hbudname )//'_BU_RRH' ) + END IF + + if ( lbudget_sv .and. (hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then + do ji = 1, 3 + call Budget( prsvs(:, :, :, nsv_c2r2beg - 1 + ji), 12 + nsv_c2r2beg - 1 + ji, Trim( hbudname )//'_BU_RSV' ) + end do + end if + if ( lbudget_sv .and. hcloud == 'LIMA' ) then + do ji = nsv_lima_beg, nsv_lima_end + call Budget( prsvs(:, :, :, ji), 12 + ji, Trim( hbudname )//'_BU_RSV' ) + 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( prths(:, :, :) * prhodj(:, :, :), 4, Trim( hbudname )//'_BU_RTH' ) + if ( lbudget_rv) call Budget( prrs (:, :, :, 1) * prhodj(:, :, :), 6, Trim( hbudname )//'_BU_RRV' ) + if ( lbudget_rc) call Budget( prrs (:, :, :, 2) * prhodj(:, :, :), 7, Trim( hbudname )//'_BU_RRC' ) + if ( lbudget_rr) call Budget( prrs (:, :, :, 3) * prhodj(:, :, :), 8, Trim( hbudname )//'_BU_RRR' ) + if ( lbudget_ri) call Budget( prrs (:, :, :, 4) * prhodj(:, :, :), 9, Trim( hbudname )//'_BU_RRI' ) + if ( lbudget_rs) call Budget( prrs (:, :, :, 5) * prhodj(:, :, :), 10, Trim( hbudname )//'_BU_RRS' ) + if ( lbudget_rg) call Budget( prrs (:, :, :, 6) * prhodj(:, :, :), 11, Trim( hbudname )//'_BU_RRG' ) + if ( lbudget_rh) call Budget( prrs (:, :, :, 7) * prhodj(:, :, :), 12, Trim( hbudname )//'_BU_RRH' ) + end if + + if ( lbudget_sv .and. (hcloud == 'C2R2' .or. hcloud == 'KHKO' )) then + do ji = 1, 3 + call Budget( prsvs(:, :, :, nsv_c2r2beg - 1 + ji), 12 + nsv_c2r2beg - 1 + ji, Trim( hbudname )//'_BU_RSV' ) + end do + end if + if ( lbudget_sv .and. hcloud == 'LIMA' ) then + do ji = nsv_lima_beg, nsv_lima_end + call Budget( prsvs(:, :, :, ji), 12 + ji, Trim( hbudname )//'_BU_RSV' ) + end do + end if +end if + +end subroutine Sources_neg_correct + +end module mode_sources_neg_correct diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90 index 500878f86..1c7aba192 100644 --- a/src/MNH/turb.f90 +++ b/src/MNH/turb.f90 @@ -342,10 +342,10 @@ END MODULE MODI_TURB !! 01/2018 (Q.Rodier) Introduction of RM17 !! 03/2020 (B.Vie) : LIMA negativity checks after turbulence, advection and microphysics budgets ! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices -! P. Wautelet + Benoit Vié 11/06/2020: improve removal of negative scalar variables + adapt the corresponding budgets -! P. Wautelet + Benoît Vié 23/06/2020: improve removal of negative scalar variables -!! -------------------------------------------------------------------------- -! +! P. Wautelet + Benoit Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets +! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct +! -------------------------------------------------------------------------- +! !* 0. DECLARATIONS ! ------------ ! @@ -380,6 +380,7 @@ USE MODI_GET_HALO USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL USE MODE_FMWRIT USE MODE_SBL +use mode_sources_neg_correct, only: Sources_neg_correct ! USE MODI_EMOIST USE MODI_ETHETA @@ -1120,179 +1121,10 @@ IF ( KRRL >= 1 ) THEN PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2) END IF END IF -! -SELECT CASE ( HCLOUD ) - CASE('KESS') - ZEXNE(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZTT(:,:,:)= PTHLT(:,:,:)*ZEXNE(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZTT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - WHERE (PRRS(:,:,:,2) < 0.) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0.0 - END WHERE - WHERE ((PRRS(:,:,:,1) <0.) .AND. (PRRS(:,:,:,2)> 0.) ) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0. - END WHERE -! - IF (LBUDGET_TH) CALL BUDGET (PRTHLS(:,:,:), 4,'NETUR_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NETUR_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NETUR_BU_RRC') - CASE('ICE3','ICE4') - ZEXNE(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZTT(:,:,:)= PTHLT(:,:,:)*ZEXNE(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZTT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZTT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - WHERE (PRRS(:,:,:,4) < 0.) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,4) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,4) = 0. - END WHERE -! -! cloud - WHERE (PRRS(:,:,:,2) < 0.) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0. - END WHERE -! -! 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) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0. - END WHERE -! ice - IF(KRR > 3) THEN - WHERE ((PRRS(:,:,:,1) < 0.).AND.(PRRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRRS(:,:,:,1),PRRS(:,:,:,4)) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + ZCOR(:,:,:) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,4) = PRRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF -! - IF (LBUDGET_TH) CALL BUDGET (PRTHLS(:,:,:), 4,'NETUR_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NETUR_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NETUR_BU_RRC') - IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4), 9,'NETUR_BU_RRI') -! - CASE('C2R2','KHKO') - ZEXNE(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZTT(:,:,:)= PTHLT(:,:,:)*ZEXNE(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZTT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZTT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) - WHERE (PRRS(:,:,:,2) < 0. .OR. PRSVS(:,:,:,NSV_C2R2BEG+1) < 0.) - PRSVS(:,:,:,NSV_C2R2BEG) = 0.0 - END WHERE - DO JSV = 2, 3 - WHERE (PRRS(:,:,:,JSV) < 0. .OR. PRSVS(:,:,:,NSV_C2R2BEG-1+JSV) < 0.) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,JSV) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,JSV) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,JSV) = 0.0 - PRSVS(:,:,:,NSV_C2R2BEG-1+JSV) = 0.0 - END WHERE - END DO - WHERE ((PRRS(:,:,:,1) <0.) .AND. (PRRS(:,:,:,2)> 0.) ) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0. - PRSVS(:,:,:,NSV_C2R2BEG+1) = 0. - END WHERE -! - IF (LBUDGET_TH) CALL BUDGET (PRTHLS(:,:,:), 4,'NETUR_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NETUR_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NETUR_BU_RRC') - IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3), 8,'NETUR_BU_RRR') - IF (LBUDGET_SV) THEN - DO JSV = 1, 3 - CALL BUDGET (PRSVS(:,:,:,NSV_C2R2BEG-1+JSV),12+NSV_C2R2BEG-1+JSV,'NETUR_BU_RSV') - END DO - END IF -! - CASE('LIMA') - ZEXNE(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZTT(:,:,:)= PTHLT(:,:,:)*ZEXNE(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZTT(:,:,:)-XTT) - ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZTT(:,:,:)-XTT) - ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1) -! Correction where rc<0 or Nc<0 - IF (LWARM) THEN - WHERE (PRRS(:,:,:,2) < XRTMIN(2)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NC) < XCTMIN(2)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - WHERE ((PRRS(:,:,:,1) <0.) .AND. (PRRS(:,:,:,2)> 0.) ) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,2) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,2) = 0. - PRSVS(:,:,:,NSV_LIMA_NC) = 0.0 - END WHERE - END IF -! Correction where rr<0 or Nr<0 - IF (LWARM .AND. LRAIN) THEN - WHERE (PRRS(:,:,:,3) < XRTMIN(3)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NR) < XCTMIN(3)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,3) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,3) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,3) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NR) = 0.0 - END WHERE - END IF -! Correction where ri<0 or Ni<0 - IF (LCOLD) THEN - WHERE (PRRS(:,:,:,4) < XRTMIN(4)/PTSTEP .OR. PRSVS(:,:,:,NSV_LIMA_NI) < XCTMIN(4)/PTSTEP) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,4) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,4) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,4) = 0.0 - PRSVS(:,:,:,NSV_LIMA_NI) = 0.0 - END WHERE - IF(KRR > 3) THEN - WHERE ((PRRS(:,:,:,1) < 0.).AND.(PRRS(:,:,:,4) > 0.)) - ZCOR(:,:,:)=MIN(-PRRS(:,:,:,1),PRRS(:,:,:,4)) - PRRS(:,:,:,1) = PRRS(:,:,:,1) + ZCOR(:,:,:) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,4) = PRRS(:,:,:,4) -ZCOR(:,:,:) - END WHERE - END IF - END IF -! - PRSVS(:, :, :, NSV_LIMA_BEG:NSV_LIMA_END) = MAX( 0.0, PRSVS(:, :, :, NSV_LIMA_BEG:NSV_LIMA_END) ) -! - IF (LBUDGET_TH) CALL BUDGET (PRTHLS(:,:,:), 4,'NETUR_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NETUR_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NETUR_BU_RRC') - IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3), 8,'NETUR_BU_RRR') - IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4), 9,'NETUR_BU_RRI') - IF (LBUDGET_SV) THEN - DO JSV = NSV_LIMA_BEG, NSV_LIMA_END - CALL BUDGET (PRSVS(:,:,:,JSV),12+JSV,'NETUR_BU_RSV') - END DO - END IF -! -END SELECT -! +! 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 ) + !---------------------------------------------------------------------------- ! !* 9. LES averaged surface fluxes -- GitLab