diff --git a/src/SURFEX/hydro.F90 b/src/SURFEX/hydro.F90
index fca65affef22d0f38e8db74af90ea1597370ac63..c77c31f80850e56167b215401219ca5851c2b5d3 100644
--- a/src/SURFEX/hydro.F90
+++ b/src/SURFEX/hydro.F90
@@ -202,9 +202,10 @@ REAL, DIMENSION(SIZE(PVEG)) :: ZWGI_EXCESS, ZF2
REAL, DIMENSION(SIZE(PEK%XWG,1),SIZE(PEK%XWG,2)) :: ZQSAT, ZQSATI, ZTI, ZPS
! For specific humidity at saturation computation (ISBA-DIF)
!
-REAL, DIMENSION(SIZE(PEK%XWG,1),SIZE(PEK%XWG,2)) :: ZWGI0
+REAL, DIMENSION(SIZE(PEK%XWG,1),SIZE(PEK%XWG,2)) :: ZWGI0, ZPHASE
! ZWGI0 = initial soil ice content (m3 m-3) before update
-! for budget diagnostics
+! for budget diagnostics
+! ZPHASE= phase change (freeze-thaw) energy (W/m2) by layer
!
!* 0.3 declarations of local parameters
!
@@ -253,7 +254,8 @@ PDELPHASEG(:) = 0.0
PDELPHASEG_SFC(:)= 0.0
ZWGI0(:,:) = 0.0
!
-ZF2(:) = MAX(XDENOM_MIN,PF2(:))
+ZF2(:) = MAX(XDENOM_MIN,PF2(:))
+ZPHASE(:,:) = 0.0
!
! Initialize evaporation components: variable definitions
! depend on snow or explicit canopy scheme:
@@ -490,25 +492,24 @@ IF (IO%CISBA=='DIF') THEN
PF2WGHT, PPS, ZQSAT, ZQSATI, ZDRAIN, ZHORTON, INL, ZQSB )
!
- CALL ICE_SOILDIF(KK, PK, PEK, ZTSTEP, ZKSFC_IVEG, ZLEGI, PSOILHCAPZ, ZWGI_EXCESS )
+ CALL ICE_SOILDIF(KK, PK, PEK, ZTSTEP, ZKSFC_IVEG, ZLEGI, PSOILHCAPZ, ZWGI_EXCESS, ZPHASE )
!
DEK%XDRAIN(:) = DEK%XDRAIN(:) + (ZDRAIN(:)+ZQSB(:)+ZWGI_EXCESS(:))/REAL(INDT)
DEK%XQSB (:) = DEK%XQSB (:) + ZQSB (:)/REAL(INDT)
DEK%XHORT (:) = DEK%XHORT (:) + ZHORTON(:)/REAL(INDT)
-!
-! Output diagnostics:
-! Compute latent heating from phase change only in surface layer and total soil column,
-! then adjust surface and total soil heat content to maintain balance.
-!
- PDELPHASEG_SFC(:) = (PEK%XWGI(:,1)-ZWGI0(:,1))*(XLMTT*XRHOLW/PTSTEP)*PK%XDZG(:,1) + &
- ZLEGI(:)*(XRHOLW*XLSTT)
- PDELPHASEG(:) = PDELPHASEG_SFC(:)
- DO JL=2,INL
- DO JJ=1,INJ
- PDELPHASEG(JJ) = PDELPHASEG(JJ) + (PEK%XWGI(JJ,JL)-ZWGI0(JJ,JL))*&
- (XLMTT*XRHOLW/PTSTEP)*PK%XDZG(JJ,JL)
+
+! Compute latent heating from phase change only in surface layer and total soil column
+
+ PDELPHASEG_SFC(:) = PDELPHASEG_SFC(:) + ZPHASE(:,1)/REAL(INDT)
+ DO JL=1,INL
+ DO JJ=1,INJ
+ PDELPHASEG(JJ) = PDELPHASEG(JJ) + ZPHASE(JJ,JL)/REAL(INDT)
ENDDO
ENDDO
+
+! Output diagnostics:
+! Adjust surface and total soil heat content to maintain balance.
+
PDELHEATG_SFC(:) = PDELHEATG_SFC(:) + PDELPHASEG_SFC(:)
PDELHEATG(:) = PDELHEATG(:) + PDELPHASEG(:)
diff --git a/src/SURFEX/ice_soildif.F90 b/src/SURFEX/ice_soildif.F90
index 2d3e2bdef90fcc10cce204ed085926cda6800c3c..9e6e29b445cbc0c12bd32da03549a57e8e007bfd 100644
--- a/src/SURFEX/ice_soildif.F90
+++ b/src/SURFEX/ice_soildif.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
- SUBROUTINE ICE_SOILDIF(KK, PK, PEK, PTSTEP, PKSFC_IVEG, PLEGI, PSOILHCAPZ, PWGI_EXCESS)
+ SUBROUTINE ICE_SOILDIF(KK, PK, PEK, PTSTEP, PKSFC_IVEG, PLEGI, PSOILHCAPZ, PWGI_EXCESS, PPHASE)
! ##########################################################################
!
!!**** *ICE_SOILDIF*
@@ -103,6 +103,9 @@ REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ
REAL, DIMENSION(:), INTENT(OUT) :: PWGI_EXCESS
! PWGI_EXCESS = Soil ice excess water content
!
+REAL, DIMENSION(:,:), INTENT(OUT) :: PPHASE
+! PPHASE = total phase change energy (just a diagnostic) (J/m3)
+!
!* 0.2 declarations of local variables
!
INTEGER :: JJ, JL ! loop control
@@ -116,7 +119,7 @@ REAL, DIMENSION(SIZE(PEK%XTG,1),SIZE(PEK%XTG,2)) :: ZK, ZEXCESSFC
REAL, DIMENSION(SIZE(PEK%XTG,1)) :: ZEXCESS
!
REAL :: ZWGMAX, ZPSIMAX, ZPSI, ZDELTAT, &
- ZPHASE, ZTGM, ZWGM, ZWGIM, ZLOG, &
+ ZTGM, ZWGM, ZWGIM, ZLOG, &
ZEFFIC, ZPHASEM, ZPHASEF, ZWORK, &
ZAPPHEATCAP
!
@@ -146,6 +149,7 @@ ZK(:,1) = PKSFC_IVEG(:)
! 2. Soil ice evolution computation:
! -------------------------------
!
+PPHASE(:,:) = 0.
DO JL=1,INL
DO JJ=1,INI
IDEPTH=PK%NWG_LAYER(JJ)
@@ -200,12 +204,14 @@ DO JL=1,INL
PEK%XTG(JJ,JL) = ZTGM + (ZPHASEF - ZPHASEM)/(PSOILHCAPZ(JJ,JL)+ZAPPHEATCAP)
!
! Get estimate of actual total phase change (J/m3) for equivalent soil water changes:
- ZPHASE = (PSOILHCAPZ(JJ,JL)+ZAPPHEATCAP)*(PEK%XTG(JJ,JL)-ZTGM)
+ PPHASE(JJ,JL) = (PSOILHCAPZ(JJ,JL)+ZAPPHEATCAP)*(PEK%XTG(JJ,JL)-ZTGM)
!
! Adjust ice and liquid water conents (m3/m3) accordingly :
- PEK%XWGI(JJ,JL) = ZWGIM + ZPHASE/(XLMTT*XRHOLW)
- PEK%XWG(JJ,JL) = ZWGM - ZPHASE/(XLMTT*XRHOLW)
+ PEK%XWGI(JJ,JL) = ZWGIM + PPHASE(JJ,JL)/(XLMTT*XRHOLW)
+ PEK%XWG(JJ,JL) = ZWGM - PPHASE(JJ,JL)/(XLMTT*XRHOLW)
!
+! Final unit conversion (W m-2):
+ PPHASE(JJ,JL) = PPHASE(JJ,JL)*PK%XDZG(JJ,JL)/PTSTEP
ENDIF
ENDDO
ENDDO
@@ -257,15 +263,17 @@ ENDDO
! and conserve energy:
!
DO JL=1,INL
- DO JJ=1,INI
- IDEPTH=PK%NWG_LAYER(JJ)
- IF(JL<=IDEPTH.AND.PEK%XWGI(JJ,JL)>0.0.AND.PEK%XWGI(JJ,JL)<1.0E-6)THEN
- PEK%XWG (JJ,JL) = PEK%XWG(JJ,JL) + PEK%XWGI(JJ,JL)
- ZEXCESSFC(JJ,JL) = ZEXCESSFC(JJ,JL) + PEK%XWGI(JJ,JL)
- PEK%XWGI(JJ,JL) = 0.0
- ENDIF
- PEK%XTG(JJ,JL) = PEK%XTG(JJ,JL) - ZEXCESSFC(JJ,JL)*XLMTT*XRHOLW/PSOILHCAPZ(JJ,JL)
- ENDDO
+ DO JJ=1,INI
+ IDEPTH=PK%NWG_LAYER(JJ)
+ IF(JL<=IDEPTH.AND.PEK%XWGI(JJ,JL)>0.0.AND.PEK%XWGI(JJ,JL)<1.0E-6)THEN
+ ZWGIM = PEK%XWGI(JJ,JL)
+ PEK%XWG(JJ,JL) = PEK%XWG(JJ,JL) + ZWGIM
+ PEK%XWGI(JJ,JL) = 0.0
+ ZPHASEM = ZWGIM*XLMTT*XRHOLW ! J m-3
+ PEK%XTG(JJ,JL) = PEK%XTG(JJ,JL) - ZPHASEM/PSOILHCAPZ(JJ,JL) ! K
+ PPHASE(JJ,JL) = PPHASE(JJ,JL) - ZPHASEM*PK%XDZG(JJ,JL)/PTSTEP ! W m-2
+ ENDIF
+ ENDDO
ENDDO
!
!
diff --git a/src/SURFEX/mode_read_surf_cov.F90 b/src/SURFEX/mode_read_surf_cov.F90
index 1a745624a059d80653d6c67a443f2f8b17c5be26..4948ea7564ae42b60dcb31d421eab1145cec2c02 100644
--- a/src/SURFEX/mode_read_surf_cov.F90
+++ b/src/SURFEX/mode_read_surf_cov.F90
@@ -258,7 +258,7 @@ ELSE
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (LHOOK) CALL DR_HOOK('READ_SURF_COV_4',0,ZHOOK_HANDLE_OMP)
#ifdef SFX_MPI
-!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JPROC,IDX,ISTATUS,INFOMPI) FIRSTPRIVATE(ZWORKR)
+!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JPROC,IDX,ZWORKR,IVAL,ISTATUS,INFOMPI)
#endif
DO JPROC=0,NPROC-1
!
diff --git a/src/SURFEX/mode_read_surf_layers.F90 b/src/SURFEX/mode_read_surf_layers.F90
index c1f6d78d65faf81cb1136a7f9a8c6ae6a428b121..cd5f5afc45f8a4bbe18f06f56d9ac136a8a04d70 100644
--- a/src/SURFEX/mode_read_surf_layers.F90
+++ b/src/SURFEX/mode_read_surf_layers.F90
@@ -248,7 +248,7 @@ ELSE
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (LHOOK) CALL DR_HOOK('READ_SURF_LAYERS_4',0,ZHOOK_HANDLE_OMP)
#ifdef SFX_MPI
-!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JPROC,IDX,ISTATUS,INFOMPI)
+!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JPROC,IDX,ZWORKR,IVAL,ISTATUS,INFOMPI)
#endif
DO JPROC=0,NPROC-1
!
diff --git a/src/SURFEX/prep_hor_ocean_field.F90 b/src/SURFEX/prep_hor_ocean_field.F90
index 7f654d671bfc21b431b64356b4d959db5725768a..9bb9915bdcd0024ff43f0c185d13793cab996479 100644
--- a/src/SURFEX/prep_hor_ocean_field.F90
+++ b/src/SURFEX/prep_hor_ocean_field.F90
@@ -84,7 +84,7 @@ REAL, POINTER, DIMENSION(:,:,:) ::ZFIELDIN=>NULL()!field to interpolate horiz
REAL, POINTER, DIMENSION(:,:) ::ZFIELD=>NULL() !field to interpolate horizontally
REAL, ALLOCATABLE, DIMENSION(:,:,:)::ZFIELDOUT!field interpolated horizontally
!
-INTEGER :: JLEV ! loop on oceanic vertical level
+INTEGER :: JLEV, JLEV2 ! loop on oceanic vertical level
INTEGER :: IK1
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!----------------------------------------------------------------------------
@@ -114,10 +114,11 @@ END IF
ALLOCATE(ZFIELDOUT (KLAT,SIZE(ZFIELDIN,2),SIZE(ZFIELDIN,3)) )
ALLOCATE(ZFIELD(SIZE(ZFIELDIN,1),SIZE(ZFIELDIN,3)))
!
-DO JLEV=1,SIZE(ZFIELDIN,2)
+DO JLEV=NOCKMIN,NOCKMAX
+ JLEV2 = JLEV - NOCKMIN + 1
WHERE (PSEABATHY(:)-XZHOC(JLEV)>0.) LINTERP(:) = .FALSE.
- ZFIELD(:,:)=ZFIELDIN(:,JLEV,:)
- CALL HOR_INTERPOL(DTCO, U, GCP, KLUOUT,ZFIELD,ZFIELDOUT(:,JLEV,:))
+ ZFIELD(:,:)=ZFIELDIN(:,JLEV2,:)
+ CALL HOR_INTERPOL(DTCO, U, GCP, KLUOUT,ZFIELD,ZFIELDOUT(:,JLEV2,:))
LINTERP(:) = .TRUE.
ENDDO
!
diff --git a/src/SURFEX/snow3L_isba.F90 b/src/SURFEX/snow3L_isba.F90
index 0f5173bed70874fc2ce4d610fbe78e93e595ab0c..de2ff87a39f01fd278480e89f07d8257dd2804ea 100644
--- a/src/SURFEX/snow3L_isba.F90
+++ b/src/SURFEX/snow3L_isba.F90
@@ -514,8 +514,6 @@ IF (PEK%TSNOW%SCHEME=='3-L' .OR. IO%CISBA == 'DIF' .OR. PEK%TSNOW%SCHEME == 'CRO
- DMK%XSNOWHMASS(:)/PTSTEP
ZGRNDFLUXN(:) = (ZSNOWH(:)+DMK%XSNOWHMASS(:))/PTSTEP + DMK%XGFLUXSNOW(:)
ZWORK(:) = PTSTEP * ZPSN(:) * (ZGRNDFLUXN(:) - PGRNDFLUX(:) - PFLSN_COR(:))
- PTG(:,1) = PTG(:,1) + ZWORK(:)*(1.-ZWGHT(:))*PCT(:)
- PTG(:,2) = PTG(:,2) + ZWORK(:)* ZWGHT(:) *ZC2(:)
ZWORK(:) = ZWORK(:) / PTSTEP
PDELHEATG(:) = PDELHEATG(:) + ZWORK(:)
PDELHEATG_SFC(:) = PDELHEATG_SFC(:) + ZWORK(:)