Gaellle 22/6/2018 : bug surfex 8.1

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(:)
 

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
 !
 !

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
         !

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
         !

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
 !

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(:)