From 320c4ff4bf8fd6dabb50501dd4d2c10fe77a1471 Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Mon, 28 May 2018 17:01:09 +0200
Subject: [PATCH] Gaelle 28/5/18 : bug surfex 8.1

---
 src/SURFEX/convert_cover_isba.F90        |  2 +-
 src/SURFEX/get_vegn.F90                  | 19 +++++----
 src/SURFEX/horibl_surf_init.F90          | 49 +++++++++++++++++-------
 src/SURFEX/init_veg_pgdn.F90             | 19 +++++++--
 src/SURFEX/mode_read_extern.F90          | 30 +--------------
 src/SURFEX/mode_read_netcdf_mercator.F90 |  2 +-
 src/SURFEX/mode_snow3l.F90               | 28 +++++++-------
 src/SURFEX/pgd_isba_par.F90              | 47 +++++++++++------------
 src/SURFEX/prep_grib_grid.F90            |  4 +-
 src/SURFEX/prep_grid_extern.F90          |  4 +-
 src/SURFEX/prep_hor_isba_field.F90       |  4 +-
 src/SURFEX/prep_hor_ocean_field.F90      |  9 +++--
 src/SURFEX/prep_hor_ocean_fields.F90     | 26 ++++++++++---
 src/SURFEX/read_gridtype_ign.F90         |  4 +-
 src/SURFEX/read_nam_grid_ign.F90         |  4 +-
 src/SURFEX/read_oceann.F90               | 37 +++++++++++-------
 src/SURFEX/write_field_2d_patch.F90      |  6 +--
 17 files changed, 167 insertions(+), 127 deletions(-)

diff --git a/src/SURFEX/convert_cover_isba.F90 b/src/SURFEX/convert_cover_isba.F90
index b6e9dd6f3..c7b569cf1 100644
--- a/src/SURFEX/convert_cover_isba.F90
+++ b/src/SURFEX/convert_cover_isba.F90
@@ -882,7 +882,7 @@ IF(HISBA/='DIF')THEN
   ENDDO
   !
   DO JLAYER=1,KGROUND
-     CALL AV_PGD(YDTCO, PDG(:,JLAYER,:),PCOVER,ZDATA_DG(:,JLAYER,:),YNAT,'ARI',OCOVER,KDECADE=KDECADE)
+     CALL AV_PGD(YDTCO, PDG(:,JLAYER,:),PCOVER,ZDATA_DG(:,JLAYER,:),YNAT,CDGAVG,OCOVER,KDECADE=KDECADE)
   ENDDO
   !
 ELSE
diff --git a/src/SURFEX/get_vegn.F90 b/src/SURFEX/get_vegn.F90
index 7a0237ff8..a41238194 100644
--- a/src/SURFEX/get_vegn.F90
+++ b/src/SURFEX/get_vegn.F90
@@ -73,9 +73,9 @@ REAL, DIMENSION(KI), INTENT(OUT) :: PLAI
 !
 TYPE(ISBA_P_t), POINTER :: PK
 TYPE(ISBA_PE_t), POINTER :: PEK
-INTEGER                               :: JI,JJ           ! loop index over tiles
+INTEGER                               :: JI, JJ           ! loop index over tiles
 INTEGER                               :: ILUOUT       ! unit numberi
-REAL, DIMENSION(U%NSIZE_NATURE)    :: ZH_TREE, ZLAI,ZWORK
+REAL, DIMENSION(U%NSIZE_NATURE)    :: ZH_TREE, ZLAI, ZWORK
 INTEGER:: IPATCH_TRBE, IPATCH_TRBD, IPATCH_TEBE, IPATCH_TEBD, IPATCH_TENE, &
           IPATCH_BOBD, IPATCH_BONE, IPATCH_BOND, IMASK, JP
 ! 
@@ -100,12 +100,13 @@ IPATCH_BONE = VEGTYPE_TO_PATCH(NVT_BONE, IO%NPATCH)
 IPATCH_BOND = VEGTYPE_TO_PATCH(NVT_BOND, IO%NPATCH)
 
 
-ZWORK(:) = S%XVEGTYPE(:,NVT_TRBE) + S%XVEGTYPE(:,NVT_TRBD) + S%XVEGTYPE(:,NVT_TEBE) + &
-           S%XVEGTYPE(:,NVT_TEBD) + S%XVEGTYPE(:,NVT_TENE) + S%XVEGTYPE(:,NVT_BOBD) + &
-           S%XVEGTYPE(:,NVT_BONE) + S%XVEGTYPE(:,NVT_BOND)
+!ZWORK(:) = S%XVEGTYPE(:,NVT_TRBE) + S%XVEGTYPE(:,NVT_TRBD) + S%XVEGTYPE(:,NVT_TEBE) + &
+!           S%XVEGTYPE(:,NVT_TEBD) + S%XVEGTYPE(:,NVT_TENE) + S%XVEGTYPE(:,NVT_BOBD) + &
+!           S%XVEGTYPE(:,NVT_BONE) + S%XVEGTYPE(:,NVT_BOND)
 
 ZH_TREE(:) = 0.
 ZLAI(:) = 0.
+ZWORK(:) = 0.
 !
 DO JP = 1,IO%NPATCH
   !
@@ -119,12 +120,14 @@ DO JP = 1,IO%NPATCH
       !
       IMASK = PK%NR_P(JJ)
       !
-      IF (S%XVEGTYPE(IMASK,JP)/=0) THEN
+      IF (PK%XH_TREE(JJ)/=XUNDEF) THEN
         !
         ZH_TREE(IMASK) = ZH_TREE(IMASK) + PK%XH_TREE(JJ) * PK%XPATCH(JJ)
         !
         ZLAI(IMASK)  = ZLAI(IMASK) + PEK%XLAI(JJ) * PK%XPATCH(JJ)
         !
+        ZWORK(IMASK) = ZWORK(IMASK) + PK%XPATCH(JJ)
+        !
       ENDIF
       !
     ENDDO
@@ -138,7 +141,9 @@ WHERE(ZWORK(:)/=0.)
   ZLAI(:) = ZLAI(:)/ZWORK(:)
 END WHERE
 !
-ZLAI(:) = U%XNATURE(:) * ZLAI(:)
+!DO JJ = 1,SIZE(ZLAI)
+!  ZLAI(JJ) = U%XNATURE(U%NR_NATURE(JJ)) * ZLAI(JJ)
+!ENDDO
 !
 !*       2. Envoi les variables vers mesonH 
 !             ------------------------------
diff --git a/src/SURFEX/horibl_surf_init.F90 b/src/SURFEX/horibl_surf_init.F90
index 5173bc4f1..caacee50c 100644
--- a/src/SURFEX/horibl_surf_init.F90
+++ b/src/SURFEX/horibl_surf_init.F90
@@ -4,7 +4,7 @@
 !SFX_LIC for details. version 1.
 !     #########
     SUBROUTINE HORIBL_SURF_INIT(PILA1,PILO1,PILA2,PILO2,KINLA,KINLO,KOLEN,&
-                           PXOUT,PYOUT,OINTERP,OGLOBLON,OGLOBN,OGLOBS,&
+                           PXOUT,PYOUT,OINTERP,OGAUSS,OGLOBLON,OGLOBN,OGLOBS,&
                            KO,KINLO_OUT,POLA,POLO,PILO1_OUT,&
                            PILO2_OUT,PLA,PILATARRAY )  
 !   ###########################################################################
@@ -121,6 +121,7 @@ INTEGER,                   INTENT(IN)  :: KOLEN   ! size of output array
 REAL,    DIMENSION(:), INTENT(IN)  :: PXOUT   ! X (lon.) of output points
 REAL,    DIMENSION(:), INTENT(IN)  :: PYOUT   ! Y (lat.) of output points
 LOGICAL, DIMENSION(:), INTENT(IN)  :: OINTERP ! .true. where physical value is needed
+LOGICAL, INTENT(IN) :: OGAUSS
 ! 
 LOGICAL, INTENT(OUT)  :: OGLOBLON  ! True if the map is circular
 LOGICAL, INTENT(OUT)  :: OGLOBN    ! True if the map has the north pole
@@ -143,12 +144,14 @@ REAL, DIMENSION(:), ALLOCATABLE    :: ZIDLAT   ! Deltai latitude
 REAL                               :: ZIDLA    ! Delta latitude
 REAL                               :: ZSOUTHPOLE! south pole latitude (-90 or  90)
 REAL                               :: ZNORTHPOLE! north pole latitude ( 90 or -90)
+REAL                               :: ZB1, ZB2
 !
 ! Variables implied in the extension procedure
 !
 INTEGER                            :: IOFFSET   ! Offset in map
 INTEGER                            :: IINLA     ! Number of parallel
  ! Loop counters
+INTEGER                            :: JLAT0, IOS_SAVE
 INTEGER                            :: JOPOS     ! Output position
 INTEGER                            :: JL, JL2   ! Dummy counter
 !
@@ -260,6 +263,7 @@ END IF
 POLA(:) = 0.
 POLO(:) = 0.
 !
+IOS_SAVE = 1
 DO JL = 1, KOLEN
   !
   IF (.NOT. OINTERP(JL)) CYCLE
@@ -273,17 +277,36 @@ DO JL = 1, KOLEN
   ! 3.1.1. find positions of latitudes
   IF (PRESENT(PILATARRAY)) THEN
     !
-    DO JL2 = 1,KINLA
-      IF((POLA(JL)>=PILATARRAY(JL2)-ZIDLAT(JL2)/2..AND.POLA(JL)<PILATARRAY(JL2)+ZIDLAT(JL2+1)/2.).OR.&
-         (POLA(JL)<=PILATARRAY(JL2)-ZIDLAT(JL2)/2..AND.POLA(JL)>PILATARRAY(JL2)+ZIDLAT(JL2+1)/2.)) THEN
-        KO(JL,3) = JL2
-        EXIT
-      ELSEIF (POLA(JL)>MAXVAL(PILATARRAY(:))) THEN
-        KO(JL,3) = MAXLOC(PILATARRAY,1)
-      ELSEIF (POLA(JL)<MINVAL(PILATARRAY(:))) THEN
-        KO(JL,3) = MINLOC(PILATARRAY,1)
-      ENDIF
-    ENDDO
+    IF (POLA(JL)>MAXVAL(PILATARRAY(:))) THEN
+      KO(JL,3) = MAXLOC(PILATARRAY,1)
+    ELSEIF (POLA(JL)<MINVAL(PILATARRAY(:))) THEN
+      KO(JL,3) = MINLOC(PILATARRAY,1)
+    ELSE 
+      DO JLAT0 = 1,KINLA
+        JL2 = MOD(JLAT0+IOS_SAVE-1,KINLA)
+        IF (JL2==0) JL2 = KINLA
+        IF (OGAUSS) THEN
+          IF (POLA(JL)>=PILATARRAY(JL2).AND.POLA(JL)<PILATARRAY(JL2+1)) THEN
+            KO(JL,3) = JL2
+            EXIT
+          ELSEIF (POLA(JL)<=PILATARRAY(JL2).AND.POLA(JL)>PILATARRAY(JL2+1)) THEN
+            KO(JL,3) = JL2
+            EXIT
+          ENDIF                 
+        ELSE
+          IF (POLA(JL)>=PILATARRAY(JL2)-ZIDLAT(JL2)/2..AND.POLA(JL)<PILATARRAY(JL2)+ZIDLAT(JL2+1)/2.) THEN
+            KO(JL,3) = JL2
+            EXIT
+          ELSEIF (POLA(JL)<=PILATARRAY(JL2)-ZIDLAT(JL2)/2..AND.POLA(JL)>PILATARRAY(JL2)+ZIDLAT(JL2+1)/2.) THEN
+            KO(JL,3) = JL2
+            EXIT
+          ENDIF                  
+        ENDIF
+              
+      ENDDO
+    ENDIF
+    IOS_SAVE = KO(JL,3)
+    !
     PLA(JL,3) = PILATARRAY(KO(JL,3))
     !
   ELSE
@@ -309,7 +332,7 @@ DO JL = 1, KOLEN
     ELSE
       PLA(JL,2) = PILATARRAY(KO(JL,3)+1)
     ENDIF
-    IF (KO(JL,3)>=KINLA) THEN
+    IF (KO(JL,3)>=KINLA-1) THEN
       PLA(JL,1) = PILATARRAY(KINLA) + 2.*ZIDLAT(KINLA)
     ELSE
       PLA(JL,1) = PILATARRAY(KO(JL,3)+2)
diff --git a/src/SURFEX/init_veg_pgdn.F90 b/src/SURFEX/init_veg_pgdn.F90
index 9a12d07d1..b21e81e09 100644
--- a/src/SURFEX/init_veg_pgdn.F90
+++ b/src/SURFEX/init_veg_pgdn.F90
@@ -54,7 +54,7 @@ USE MODD_AGRI_n, ONLY : AGRI_t
 USE MODD_SURF_ATM,       ONLY : LCPL_ARP
 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
 USE MODD_SURF_PAR,       ONLY : XUNDEF, NUNDEF
-USE MODD_CSTS,           ONLY : XCPD, XLVTT, XLSTT
+USE MODD_CSTS,           ONLY : XCPD, XLVTT, XLSTT, XDAY
 USE MODD_SNOW_PAR,       ONLY : XEMISSN
 USE MODD_ISBA_PAR,       ONLY : XTAU_ICE
 !
@@ -122,9 +122,11 @@ REAL, DIMENSION(:), INTENT(IN) :: PRHOA
 !*       0.2   Declarations of local variables
 !              -------------------------------
 !
+REAL, DIMENSION(KI,IO%NGROUND_LAYER) :: ZCONDSAT
+!
 INTEGER :: JPATCH  ! loop counter on tiles
 INTEGER :: JILU,JP, JMAXLOC    ! loop increment
-INTEGER :: JL  ! loop counter on layers
+INTEGER :: JL, JI  ! loop counter on layers
 !
 INTEGER :: IABC
 !
@@ -216,7 +218,17 @@ IF (.NOT.ASSOCIATED(K%XMPOTSAT)) THEN
     ALLOCATE(K%XWD0   (KI,IO%NGROUND_LAYER))
     ALLOCATE(K%XKANISO(KI,IO%NGROUND_LAYER))
     !
-    IF(IO%CISBA=='DIF')THEN
+    IF (DTI%LDATA_WFC.OR.DTI%LDATA_WSAT) THEN
+      IF (DTI%LDATA_CONDSAT) THEN
+        ZCONDSAT(:,:) = DTI%XPAR_CONDSAT(:,:)
+      ELSE
+        DO JL = 1,IO%NGROUND_LAYER
+          ZCONDSAT(:,JL) = HYDCONDSAT_FUNC(K%XCLAY(:,JL),K%XSAND(:,JL),IO%CPEDOTF)
+        ENDDO
+      ENDIF
+      K%XWD0(:,:) = K%XWSAT(:,:) * ((0.0001/XDAY)/ZCONDSAT(:,:))**(1./(2.*K%XBCOEF(:,:)+3.))
+      print*,'wd0 ',minval(K%XWD0),maxval(K%XWD0)
+    ELSEIF(IO%CISBA=='DIF')THEN
       K%XWD0(:,:) = WFC_FUNC(K%XCLAY(:,:),K%XSAND(:,:),IO%CPEDOTF)
     ELSE
       K%XWD0(:,:) = K%XWWILT(:,:)
@@ -469,6 +481,7 @@ ELSE
     PK%XCONDSAT(:,JL) = HYDCONDSAT_FUNC(KK%XCLAY(:,JL),KK%XSAND(:,JL),IO%CPEDOTF) 
   END DO
 ENDIF
+!
 PK%XTAUICE(:) = XTAU_ICE
 !
 IF (IO%CISBA=='2-L' .OR. IO%CISBA=='3-L') THEN
diff --git a/src/SURFEX/mode_read_extern.F90 b/src/SURFEX/mode_read_extern.F90
index 96f488f92..a5bbd8586 100644
--- a/src/SURFEX/mode_read_extern.F90
+++ b/src/SURFEX/mode_read_extern.F90
@@ -644,7 +644,7 @@ IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ') THEN
     ALLOCATE(ZCONDSAT(KNI,ILAYER))    
   ENDIF 
   !
-  IF (GDIM) THEN
+  IF (GDIM.AND..NOT.GTEB) THEN
     YRECFM='L_WFC'
     CALL READ_SURF(HFILEPGDTYPE,YRECFM,GDATA_WFC,IRESP,HDIR='E')
     IF (GDATA_WFC) THEN
@@ -668,33 +668,7 @@ IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ') THEN
         WRITE(YRECFM,FMT='(A9,I2.2)') 'D_WSAT_L',JL
         CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZWSAT(:,JL),IRESP,HDIR='E')
       ENDDO 
-    ENDIF
-    IF (GTEB.AND..NOT.GGD) THEN         
-      YRECFM='L_CONDSAT'
-      CALL READ_SURF(HFILEPGDTYPE,YRECFM,GDATA_CONDSAT,IRESP,HDIR='E')
-      IF (GDATA_CONDSAT) THEN
-        DO JL=1,ILAYER
-          WRITE(YRECFM,FMT='(A9,I2.2)') 'D_CNDSAT_L',JL
-          CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZCONDSAT(:,JL),IRESP,HDIR='E')
-        ENDDO 
-      ENDIF
-      YRECFM='L_MPOTSAT'
-      CALL READ_SURF(HFILEPGDTYPE,YRECFM,GDATA_MPOTSAT,IRESP,HDIR='E')
-      IF (GDATA_MPOTSAT) THEN
-        DO JL=1,ILAYER
-          WRITE(YRECFM,FMT='(A9,I2.2)') 'D_MPTSAT_L',JL
-          CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZMPOTSAT(:,JL),IRESP,HDIR='E')
-        ENDDO 
-      ENDIF
-      YRECFM='L_BCOEF'
-      CALL READ_SURF(HFILEPGDTYPE,YRECFM,GDATA_BCOEF,IRESP,HDIR='E')
-      IF (GDATA_BCOEF) THEN
-        DO JL=1,ILAYER
-          WRITE(YRECFM,FMT='(A9,I2.2)') 'D_BCOEF_L',JL
-          CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZBCOEF(:,JL),IRESP,HDIR='E')
-        ENDDO 
-      ENDIF
-    ENDIF      
+    ENDIF    
   ELSE
     GDATA_WWILT   = .FALSE.
     GDATA_WFC     = .FALSE.
diff --git a/src/SURFEX/mode_read_netcdf_mercator.F90 b/src/SURFEX/mode_read_netcdf_mercator.F90
index b1469fd3b..7af093418 100644
--- a/src/SURFEX/mode_read_netcdf_mercator.F90
+++ b/src/SURFEX/mode_read_netcdf_mercator.F90
@@ -771,7 +771,7 @@ IF (ALLOCATED(XLAT_OUT)) THEN
   IINLA = NINLAT
   ALLOCATE(NINLOH(IINLA+4))
   CALL HORIBL_SURF_INIT(XILAT1,XILON1,XILAT2,XILON2,NINLAT,NINLON, &
-                        INO,XLON_OUT,XLAT_OUT,LINTERP,LGLOBLON,&
+                        INO,XLON_OUT,XLAT_OUT,LINTERP,.FALSE.,LGLOBLON,&
                         LGLOBN,LGLOBS,NO,NINLOH,XOLA,XOLO,XILO1H,&
                         XILO2H,XLA,XILATARRAY)
   !
diff --git a/src/SURFEX/mode_snow3l.F90 b/src/SURFEX/mode_snow3l.F90
index 85a55f9a6..c62b91ba4 100644
--- a/src/SURFEX/mode_snow3l.F90
+++ b/src/SURFEX/mode_snow3l.F90
@@ -1165,16 +1165,16 @@ ELSE !(INLVLS>=10 and /=12)
 !
 ENDIF
 !
-DO JJ=1,INLVLS
-   DO JI=1,INI
-      IF(PSNOW(JI)==XUNDEF)THEN
-         PSNOWDZ(JI,JJ) = XUNDEF
-      ELSEIF(.NOT.GREGRID(JI))THEN           
-         PSNOWDZ(JI,JJ)=PSNOWDZ_OLD(JI,JJ)
-      ENDIF      
-   ENDDO
+DO JI=1,INI
+  IF(PSNOW(JI)==XUNDEF) PSNOWDZ(JI,:) = XUNDEF
 ENDDO
 !
+IF (PRESENT(PSNOWDZ_OLD)) THEN
+  DO JI=1,INI
+    IF (.NOT.GREGRID(JI)) PSNOWDZ(JI,:)=PSNOWDZ_OLD(JI,:) 
+  ENDDO
+ENDIF
+!
 IF (LHOOK) CALL DR_HOOK('MODE_SNOW3L:SNOW3LGRID_2D',1,ZHOOK_HANDLE)
 !
 END SUBROUTINE SNOW3LGRID_2D
@@ -1463,13 +1463,11 @@ ELSE
 !
 ENDIF
 !
-DO JJ=1,INLVLS
-  IF(PSNOW==XUNDEF)THEN
-    PSNOWDZ(JJ) = XUNDEF
-  ELSEIF(.NOT.GREGRID)THEN
-    PSNOWDZ(JJ) = PSNOWDZ_OLD(JJ)
-  ENDIF
-END DO
+IF (PSNOW==XUNDEF) PSNOWDZ(:) = XUNDEF
+!
+IF (PRESENT(PSNOWDZ_OLD)) THEN
+  IF (.NOT.GREGRID) PSNOWDZ(:) = PSNOWDZ_OLD(:)
+ENDIF
 !
 IF (LHOOK) CALL DR_HOOK('MODE_SNOW3L:SNOW3LGRID_1D',1,ZHOOK_HANDLE)
 !
diff --git a/src/SURFEX/pgd_isba_par.F90 b/src/SURFEX/pgd_isba_par.F90
index 3bb5202c2..9c8ac6ff6 100644
--- a/src/SURFEX/pgd_isba_par.F90
+++ b/src/SURFEX/pgd_isba_par.F90
@@ -1154,7 +1154,28 @@ IF (IO%CPHOTO/='NON') THEN
                              HPROGRAM,'MAJ','REAP_D: day of reaping','NAT',&
                              CFNAM_REAP_D,CFTYP_REAP_D,XUNIF_REAP_D,DTV%XPAR_REAP_D,DTV%LDATA_REAP_D)
       IF (ALL(.NOT.DTV%LDATA_REAP_D)) DEALLOCATE(DTV%XPAR_REAP_D)
-      !      
+      !  
+      IF (.NOT.IO%LECOCLIMAP .AND. .NOT.(ANY(DTV%LDATA_IRRIG).AND.ANY(DTV%LDATA_WATSUP).AND.&
+            ANY(DTV%LDATA_SEED_M).AND.ANY(DTV%LDATA_SEED_D).AND.ANY(DTV%LDATA_REAP_M).AND.&
+            ANY(DTV%LDATA_REAP_D))) THEN
+        !
+        WRITE(ILUOUT,*) ' '
+        WRITE(ILUOUT,*) '***********************************************************'
+        WRITE(ILUOUT,*) '* Error in PGD field preparation of irrigation fields      *'
+        WRITE(ILUOUT,*) '* There is no prescribed value and no input file :         *'
+        IF (ALL(.NOT.DTV%LDATA_IRRIG   )) WRITE(ILUOUT,*) '* for IRRIG                   *'
+        IF (ALL(.NOT.DTV%LDATA_WATSUP  )) WRITE(ILUOUT,*) '* for WATSUP                  *'
+        IF (ALL(.NOT.DTV%LDATA_SEED_M  )) WRITE(ILUOUT,*) '* for SEED_M                  *'
+        IF (ALL(.NOT.DTV%LDATA_SEED_D  )) WRITE(ILUOUT,*) '* for SEED_D                  *'
+        IF (ALL(.NOT.DTV%LDATA_REAP_M  )) WRITE(ILUOUT,*) '* for REAP_M                  *'
+        IF (ALL(.NOT.DTV%LDATA_REAP_D  )) WRITE(ILUOUT,*) '* for REAP_D                  *'      
+        WRITE(ILUOUT,*) '* Without ECOCLIMAP, these fields must be prescribed      *'
+        WRITE(ILUOUT,*) '***********************************************************'
+        WRITE(ILUOUT,*) ' '
+        CALL ABOR1_SFX('PGD_ISBA_PAR: NO PRESCRIBED VALUE NOR INPUT FILE FOR IRRIGATION PARAMETERS')
+        !
+      ENDIF
+      
     ENDIF
     !
     IF ((ANY(DTV%LDATA_IRRIG).OR.ANY(DTV%LDATA_WATSUP).OR.ANY(DTV%LDATA_SEED_M).OR.&
@@ -1173,29 +1194,7 @@ IF (IO%CPHOTO/='NON') THEN
       WRITE(ILUOUT,*) ' '
       CALL ABOR1_SFX('PGD_ISBA_PAR: MISSING PRESCRIBED VALUE OR INPUT FILE FOR IRRIGATION PARAMETERS')
       !
-    ENDIF
-    !
-    IF (.NOT.IO%LECOCLIMAP .AND. .NOT.(ANY(DTV%LDATA_IRRIG).AND.ANY(DTV%LDATA_WATSUP).AND.&
-            ANY(DTV%LDATA_SEED_M).AND.ANY(DTV%LDATA_SEED_D).AND.ANY(DTV%LDATA_REAP_M).AND.&
-            ANY(DTV%LDATA_REAP_D))) THEN
-      !
-      WRITE(ILUOUT,*) ' '
-      WRITE(ILUOUT,*) '***********************************************************'
-      WRITE(ILUOUT,*) '* Error in PGD field preparation of irrigation fields      *'
-      WRITE(ILUOUT,*) '* There is no prescribed value and no input file :         *'
-      IF (ALL(.NOT.DTV%LDATA_IRRIG   )) WRITE(ILUOUT,*) '* for IRRIG                   *'
-      IF (ALL(.NOT.DTV%LDATA_WATSUP  )) WRITE(ILUOUT,*) '* for WATSUP                  *'
-      IF (ALL(.NOT.DTV%LDATA_SEED_M  )) WRITE(ILUOUT,*) '* for SEED_M                  *'
-      IF (ALL(.NOT.DTV%LDATA_SEED_D  )) WRITE(ILUOUT,*) '* for SEED_D                  *'
-      IF (ALL(.NOT.DTV%LDATA_REAP_M  )) WRITE(ILUOUT,*) '* for REAP_M                  *'
-      IF (ALL(.NOT.DTV%LDATA_REAP_D  )) WRITE(ILUOUT,*) '* for REAP_D                  *'      
-      WRITE(ILUOUT,*) '* Without ECOCLIMAP, these fields must be prescribed      *'
-      WRITE(ILUOUT,*) '***********************************************************'
-      WRITE(ILUOUT,*) ' '
-      CALL ABOR1_SFX('PGD_ISBA_PAR: NO PRESCRIBED VALUE NOR INPUT FILE FOR IRRIGATION PARAMETERS')
-      !
-    ENDIF
-    
+    ENDIF    
     !
   ENDIF
   !
diff --git a/src/SURFEX/prep_grib_grid.F90 b/src/SURFEX/prep_grib_grid.F90
index 288682579..238385f00 100644
--- a/src/SURFEX/prep_grib_grid.F90
+++ b/src/SURFEX/prep_grib_grid.F90
@@ -639,13 +639,13 @@ IF (ALLOCATED(XLAT_OUT)) THEN
     IINLA = NINLA     
     ALLOCATE(NINLOH(IINLA+4)) 
     CALL HORIBL_SURF_INIT(XILA1,XILO1,XILA2,XILO2,NINLA,NINLO,INO,XLON,XLAT, &
-                          LINTERP,LGLOBLON,LGLOBN,LGLOBS,NO, &
+                          LINTERP,.TRUE.,LGLOBLON,LGLOBN,LGLOBS,NO, &
                           NINLOH,XOLA,XOLO,XILO1H,XILO2H,XLA)
   ELSEIF (HGRIDTYPE=='AROME     ') THEN
     IINLA = NY
     ALLOCATE(NINLOH(IINLA+4))
     CALL HORIBL_SURF_INIT(0.,0.,XY,XX,NY,NIX,INO,XZX,XZY, &
-                          LINTERP,LGLOBLON,LGLOBN,LGLOBS,NO, &
+                          LINTERP,.FALSE.,LGLOBLON,LGLOBN,LGLOBS,NO, &
                           NINLOH,XOLA,XOLO,XILO1H,XILO2H,XLA)        
   ENDIF              
 !
diff --git a/src/SURFEX/prep_grid_extern.F90 b/src/SURFEX/prep_grid_extern.F90
index 0e01f370a..0ad52b46d 100644
--- a/src/SURFEX/prep_grid_extern.F90
+++ b/src/SURFEX/prep_grid_extern.F90
@@ -149,13 +149,13 @@ IF (ALLOCATED(XLAT_OUT)) THEN
       IINLA = NINLA
       ALLOCATE(NINLOH(IINLA+4))
       CALL HORIBL_SURF_INIT(XILA1,XILO1,XILA2,XILO2,NINLA,NINLO,INO,XLON,XLAT, &
-                            LINTERP,LGLOBLON,LGLOBN,LGLOBS,NO, &
+                            LINTERP,.TRUE.,LGLOBLON,LGLOBN,LGLOBS,NO, &
                             NINLOH,XOLA,XOLO,XILO1H,XILO2H,XLA,XILATARRAY_g)
     ELSEIF (HGRIDTYPE=='LATLON    ') THEN
       IINLA = NINLAT
       ALLOCATE(NINLOH(IINLA+4))
       CALL HORIBL_SURF_INIT(XILAT1,XILON1,XILAT2,XILON2,NINLAT,NINLON, &
-                            INO,XLON_OUT,XLAT_OUT,LINTERP,LGLOBLON,&
+                            INO,XLON_OUT,XLAT_OUT,LINTERP,.FALSE.,LGLOBLON,&
                             LGLOBN,LGLOBS,NO,NINLOH,XOLA,XOLO,XILO1H,&
                             XILO2H,XLA,XILATARRAY_l)
   
diff --git a/src/SURFEX/prep_hor_isba_field.F90 b/src/SURFEX/prep_hor_isba_field.F90
index 8e2fbee39..dfb2b2dc2 100644
--- a/src/SURFEX/prep_hor_isba_field.F90
+++ b/src/SURFEX/prep_hor_isba_field.F90
@@ -54,7 +54,9 @@ USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
 USE MODD_SSO_n, ONLY : SSO_t
 USE MODD_GRID_CONF_PROJ_n, ONLY : GRID_CONF_PROJ_t
 !
-USE MODD_PREP,     ONLY : CINGRID_TYPE, CINTERP_TYPE, XZS_LS, LINTERP, CMASK
+USE MODD_PREP,     ONLY : CINGRID_TYPE, CINTERP_TYPE, XZS_LS, &
+                          XLAT_OUT, XLON_OUT, XX_OUT, XY_OUT, &
+                          LINTERP, CMASK
 USE MODD_GRID_GRIB, ONLY : CINMODEL  
 !
 USE MODD_PREP_ISBA, ONLY : XGRID_SOIL, NGRID_LEVEL, LSNOW_IDEAL,    &
diff --git a/src/SURFEX/prep_hor_ocean_field.F90 b/src/SURFEX/prep_hor_ocean_field.F90
index be0b3cd4b..7f654d671 100644
--- a/src/SURFEX/prep_hor_ocean_field.F90
+++ b/src/SURFEX/prep_hor_ocean_field.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 PREP_HOR_OCEAN_FIELD (DTCO, UG, U, GCP, O, OR, KLAT, HPROGRAM,   &
+SUBROUTINE PREP_HOR_OCEAN_FIELD (DTCO, UG, U, GCP, O, OR, KLAT, PSEABATHY, HPROGRAM,   &
                                  HFILE,HFILETYPE,KLUOUT,OUNIF,   &
                                  HSURF,HNCVARNAME                )
 !     #######################################################
@@ -42,8 +42,8 @@ USE MODD_OCEAN_REL_n, ONLY : OCEAN_REL_t
 !
 USE MODD_CSTS,           ONLY : XTT
 USE MODD_SURF_PAR,       ONLY : XUNDEF
-USE MODD_OCEAN_GRID,   ONLY : NOCKMIN,NOCKMAX
-USE MODD_PREP,           ONLY : CINGRID_TYPE, CINTERP_TYPE
+USE MODD_OCEAN_GRID,   ONLY : NOCKMIN,NOCKMAX,XZHOC
+USE MODD_PREP,           ONLY : CINGRID_TYPE, CINTERP_TYPE, LINTERP
 !
 USE MODI_PREP_OCEAN_UNIF
 USE MODI_PREP_OCEAN_NETCDF
@@ -67,6 +67,7 @@ TYPE(GRID_CONF_PROJ_t),INTENT(INOUT) :: GCP
 TYPE(OCEAN_t), INTENT(INOUT) :: O
 TYPE(OCEAN_REL_t), INTENT(INOUT) :: OR
 INTEGER, INTENT(IN) :: KLAT
+REAL, DIMENSION(:), INTENT(IN) :: PSEABATHY
 !
  CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
  CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! file name
@@ -114,8 +115,10 @@ ALLOCATE(ZFIELDOUT  (KLAT,SIZE(ZFIELDIN,2),SIZE(ZFIELDIN,3)) )
 ALLOCATE(ZFIELD(SIZE(ZFIELDIN,1),SIZE(ZFIELDIN,3)))
 !
 DO JLEV=1,SIZE(ZFIELDIN,2)
+  WHERE (PSEABATHY(:)-XZHOC(JLEV)>0.) LINTERP(:) = .FALSE.
   ZFIELD(:,:)=ZFIELDIN(:,JLEV,:)
   CALL HOR_INTERPOL(DTCO, U, GCP, KLUOUT,ZFIELD,ZFIELDOUT(:,JLEV,:))
+  LINTERP(:) = .TRUE.
 ENDDO
 !
 !*      5.     Return to historical variable
diff --git a/src/SURFEX/prep_hor_ocean_fields.F90 b/src/SURFEX/prep_hor_ocean_fields.F90
index 93975631d..2b9515b7d 100644
--- a/src/SURFEX/prep_hor_ocean_fields.F90
+++ b/src/SURFEX/prep_hor_ocean_fields.F90
@@ -89,25 +89,25 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('PREP_HOR_OCEAN_FIELDS',0,ZHOOK_HANDLE)
 YSURF='TEMP_OC'
 YNCVARNAME='temperature'
- CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, &
+ CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, PSEABATHY, &
                            HPROGRAM,HFILE,HFILETYPE,KLUOUT,OUNIF,YSURF,YNCVARNAME)
 !---------------------------------------------------------------------------
 !
 !*      4.     Treatment of oceanic salinity
 YSURF='SALT_OC'
 YNCVARNAME='salinity'
- CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, &
+ CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, PSEABATHY, &
                            HPROGRAM,HFILE,HFILETYPE,KLUOUT,OUNIF,YSURF,YNCVARNAME)
 !---------------------------------------------------------------------------
 !
 !*      5.     Treatment of oceanic current
 YSURF='UCUR_OC'
 YNCVARNAME='u'
- CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, &
+ CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, PSEABATHY, &
                            HPROGRAM,HFILE,HFILETYPE,KLUOUT,OUNIF,YSURF,YNCVARNAME)
 YSURF='VCUR_OC'
 YNCVARNAME='v'
- CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, &
+ CALL PREP_HOR_OCEAN_FIELD(DTCO, UG, U, GCP, O, OR, KLAT, PSEABATHY, &
                            HPROGRAM,HFILE,HFILETYPE,KLUOUT,OUNIF,YSURF,YNCVARNAME)
 !---------------------------------------------------------------------------
 !
@@ -154,8 +154,22 @@ IF (IL/=0) THEN
         !
         OR%XSEAU_REL(J,JLEV)  = XUNDEF
         OR%XSEAV_REL(J,JLEV)  = XUNDEF
-        !        
-      ENDIF 
+        ! 
+      ELSEIF (O%XSEAT(J,JLEV)==XUNDEF) THEN
+        O%XSEAT(J,JLEV) = O%XSEAT(J,JLEV-1)
+        OR%XSEAT_REL(J,JLEV) = OR%XSEAT_REL(J,JLEV-1)
+      ELSEIF (O%XSEAS(J,JLEV)==XUNDEF) THEN
+        O%XSEAS(J,JLEV) = O%XSEAS(J,JLEV-1)
+        OR%XSEAS_REL(J,JLEV) = OR%XSEAS_REL(J,JLEV-1)
+      ELSEIF (O%XSEAE(J,JLEV)==XUNDEF) THEN
+        O%XSEAE(J,JLEV) = O%XSEAE(J,JLEV-1)
+      ELSEIF (O%XSEAU(J,JLEV)==XUNDEF) THEN
+        O%XSEAU(J,JLEV) = O%XSEAU(J,JLEV-1)
+        OR%XSEAU_REL(J,JLEV) = OR%XSEAU_REL(J,JLEV-1)
+      ELSEIF (O%XSEAV(J,JLEV)==XUNDEF) THEN
+        O%XSEAV(J,JLEV) = O%XSEAV(J,JLEV-1)
+        OR%XSEAV_REL(J,JLEV) = OR%XSEAV_REL(J,JLEV-1)
+      ENDIF        
     ENDDO
   ENDDO
 !
diff --git a/src/SURFEX/read_gridtype_ign.F90 b/src/SURFEX/read_gridtype_ign.F90
index c06bb243b..6b8cb4a49 100644
--- a/src/SURFEX/read_gridtype_ign.F90
+++ b/src/SURFEX/read_gridtype_ign.F90
@@ -74,8 +74,8 @@ REAL, DIMENSION(KLU)              :: ZY       ! Y  Lambert coordinate of grid me
 REAL, DIMENSION(KLU)              :: ZDX      ! X grid mesh size
 REAL, DIMENSION(KLU)              :: ZDY      ! Y grid mesh size
 !
-REAL, DIMENSION(KLU*3)            :: ZXALL    ! maximum domain X coordinate of grid mesh
-REAL, DIMENSION(KLU*3)            :: ZYALL    ! maximum domain Y coordinate of grid mesh
+REAL, DIMENSION(KLU*5)            :: ZXALL    ! maximum domain X coordinate of grid mesh
+REAL, DIMENSION(KLU*5)            :: ZYALL    ! maximum domain Y coordinate of grid mesh
 INTEGER                           :: IDIMX    ! maximum domain length in X
 INTEGER                           :: IDIMY    ! maximum domain length in Y
 INTEGER                           :: ILUOUT
diff --git a/src/SURFEX/read_nam_grid_ign.F90 b/src/SURFEX/read_nam_grid_ign.F90
index ae96b0429..f1d9ec56a 100644
--- a/src/SURFEX/read_nam_grid_ign.F90
+++ b/src/SURFEX/read_nam_grid_ign.F90
@@ -228,8 +228,8 @@ IF (HDIR/='H') THEN
   !*       7.    maximum domain lengths
   !              ----------------------
   !
-  ALLOCATE(ZXALL(KL*3))
-  ALLOCATE(ZYALL(KL*3))
+  ALLOCATE(ZXALL(KL*5))
+  ALLOCATE(ZYALL(KL*5))
   CALL GET_XYALL_IGN(ZX,ZY,ZDX,ZDY,ZXALL,ZYALL,IDIMX,IDIMY)
   !
   !---------------------------------------------------------------------------
diff --git a/src/SURFEX/read_oceann.F90 b/src/SURFEX/read_oceann.F90
index ac9e6af5f..36846ac69 100644
--- a/src/SURFEX/read_oceann.F90
+++ b/src/SURFEX/read_oceann.F90
@@ -53,7 +53,7 @@ USE MODD_OCEAN_GRID, ONLY : NOCKMIN,NOCKMAX,XZHOC
 !
 USE MODI_READ_SURF
 USE MODI_OCEAN_MERCATORVERGRID
-!
+USE MODI_PREP_OCEAN_MERCATORVERGRID
 !
 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
 USE PARKIND1  ,ONLY : JPRB
@@ -82,7 +82,7 @@ INTEGER           :: IRESP          ! Error code after redding
 !
  CHARACTER(LEN=4)  :: YLVL
 !
- CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
+ CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
  CHARACTER(LEN=14) :: YFORM          ! Writing format
 REAL, DIMENSION(:),ALLOCATABLE  :: ZWORK      ! 1D array to write data in file
 !
@@ -134,18 +134,27 @@ ENDIF
 !
 !-------------------------------------------------------------------------------
 !
-NOCKMIN = 0
-YRECFM='SEA_NBLEVEL'
- CALL READ_SURF(HPROGRAM,YRECFM,NOCKMAX,IRESP)
-!
-ALLOCATE(XZHOC(NOCKMIN:NOCKMAX))
-XZHOC(NOCKMIN) = 0.
-! Read vertical grid
-DO JLEVEL = NOCKMIN+1,NOCKMAX
-  WRITE(YLVL,'(I4)') JLEVEL
-  YRECFM='LEVL_OC'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
-  CALL READ_SURF(HPROGRAM,YRECFM,XZHOC(JLEVEL),IRESP)
-END DO
+IF (IVERSION>=8) THEN
+
+  NOCKMIN = 0
+  YRECFM='SEA_NBLEVEL'
+  CALL READ_SURF(HPROGRAM,YRECFM,NOCKMAX,IRESP)
+  !
+  ALLOCATE(XZHOC(NOCKMIN:NOCKMAX))
+  XZHOC(NOCKMIN) = 0.
+  ! Read vertical grid
+  DO JLEVEL = NOCKMIN+1,NOCKMAX
+    WRITE(YLVL,'(I4)') JLEVEL
+    YRECFM='LEVL_OC'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
+    CALL READ_SURF(HPROGRAM,YRECFM,XZHOC(JLEVEL),IRESP)
+  END DO
+  !
+ELSE
+  !
+  NOCKMIN = 0
+  CALL PREP_OCEAN_MERCATORVERGRID(HPROGRAM,.TRUE.)
+  !
+ENDIF
 !
  CALL OCEAN_MERCATORVERGRID
 !
diff --git a/src/SURFEX/write_field_2d_patch.F90 b/src/SURFEX/write_field_2d_patch.F90
index bda5f22c9..f296a11cd 100644
--- a/src/SURFEX/write_field_2d_patch.F90
+++ b/src/SURFEX/write_field_2d_patch.F90
@@ -53,11 +53,11 @@ ELSE
   !
   IF (KP/=0) THEN
     PWORK_WR(:,:,KP) = ZWORK(:,:)
-    IF ( KP==SIZE(PWORK_WR,2) ) THEN
-      CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,PWORK_WR,IRESP,HCOMMENT=HCOMMENT)
+    IF ( KP==SIZE(PWORK_WR,3) ) THEN
+      CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,PWORK_WR,IRESP,HCOMMENT=HCOMMENT,HNAM_DIM=HNAM_DIM)
    ENDIF
   ELSE
-    CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,ZWORK,IRESP,HCOMMENT=HCOMMENT)
+    CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,ZWORK,IRESP,HCOMMENT=HCOMMENT,HNAM_DIM=HNAM_DIM)
   ENDIF
   !
 ENDIF
-- 
GitLab