From a9a859fe264c9e9730fe44d24af2975034b23588 Mon Sep 17 00:00:00 2001
From: Gaelle Tanguy <gaelle.tanguy@meteo.fr>
Date: Wed, 9 Apr 2014 09:37:34 +0000
Subject: [PATCH] Gaelle 9/4/2014 : bug surfex 7.3

---
 src/SURFEX/ch_init_snapn.F90           |  1 +
 src/SURFEX/coare30_flux.F90            | 32 +++++++++------
 src/SURFEX/compute_isba_parameters.F90 | 12 +++---
 src/SURFEX/ecoclimap2_lai.F90          |  1 -
 src/SURFEX/init_isban.F90              |  1 +
 src/SURFEX/isba_snow_agr.F90           |  2 +-
 src/SURFEX/modd_surfex_omp.F90         |  4 +-
 src/SURFEX/mode_psychro.F90            |  3 +-
 src/SURFEX/mode_read_surf_asc.F90      |  2 +-
 src/SURFEX/mode_read_surf_fa.F90       | 56 ++++++++++++--------------
 src/SURFEX/mode_soil.F90               | 10 +++--
 src/SURFEX/pgd_surf_atm.F90            |  1 +
 src/SURFEX/pgd_teb_veg.F90             |  1 +
 src/SURFEX/prep_isba_netcdf.F90        |  3 +-
 src/SURFEX/prep_ocean_ascllv.F90       |  3 +-
 src/SURFEX/prep_snow_extern.F90        | 43 +++++++++++++++++---
 src/SURFEX/prep_ver_snow.F90           |  7 ++--
 src/SURFEX/prep_ver_teb.F90            |  4 +-
 src/SURFEX/read_gridtype_gauss.F90     |  6 ++-
 src/SURFEX/read_nam_grid_gauss.F90     | 12 ++++--
 src/SURFEX/snowcro.F90                 | 40 ++++++++----------
 src/SURFEX/sunpos.F90                  |  3 +-
 22 files changed, 149 insertions(+), 98 deletions(-)

diff --git a/src/SURFEX/ch_init_snapn.F90 b/src/SURFEX/ch_init_snapn.F90
index 2af97ddaf..912531b65 100644
--- a/src/SURFEX/ch_init_snapn.F90
+++ b/src/SURFEX/ch_init_snapn.F90
@@ -24,6 +24,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original        11/2011
+!!      J.Escobar       11/2013 : ajout use MODI_CH_OPEN_INPUTB
 !!-----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
diff --git a/src/SURFEX/coare30_flux.F90 b/src/SURFEX/coare30_flux.F90
index 0e853b550..50b43f62b 100644
--- a/src/SURFEX/coare30_flux.F90
+++ b/src/SURFEX/coare30_flux.F90
@@ -50,13 +50,17 @@
 !!      Original     1/06/2006
 !!      B. Decharme    06/2009 limitation of Ri
 !!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
+!!      B. Decharme    06/2013 bug in z0 (output) computation 
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
 !               ------------
 !
 USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
-                              XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT,XP00  
+                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
+                            XP00
+USE MODD_SURF_ATM,   ONLY : XVZ0CM
+!
 USE MODD_SEAFLUX_n
 USE MODD_SURF_PAR,   ONLY : XUNDEF, XSURF_EPSILON
 USE MODD_WATER_PAR
@@ -226,10 +230,12 @@ ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
 ZQASAT(:) = QSAT(PTA(:),ZPA(:)) 
 !
 !
-ZO(:)  = 0.0001
+ZO(:)  = PZ0SEA(:)
 ZWG(:) = 0.
 IF (LPWG) ZWG(:) = 0.5
 !
+ZCHARN(:) = 0.011  
+!
 DO J=1,SIZE(PTA)
   !
   !      2.2       initial guess
@@ -239,16 +245,16 @@ DO J=1,SIZE(PTA)
   ZDT(J) = -(PTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference
   ZDQ(J) = PQSAT(J)-PQA(J)                         !specific humidity difference
   !
-  ZDUWG(J) = (ZDU(J)*ZDU(J)+ZWG(J)*ZWG(J))**0.5    !wind speed difference including gustiness ZWG
+  ZDUWG(J) = SQRT(ZDU(J)**2+ZWG(J)**2)     !wind speed difference including gustiness ZWG
   !
   !      2.3   initialization of neutral coefficients
   !
   ZU10(J)  = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J))
   ZUSR(J)  = 0.035*ZU10(J)
   ZVISA(J) = 1.326E-5*(1.+6.542E-3*(PTA(J)-XTT)+&
-             8.301E-6*(PTA(J)-XTT)**2.-4.84E-9*(PTA(J)-XTT)**3.) !Andrea (1989) CRREL Rep. 89-11
+             8.301E-6*(PTA(J)-XTT)**2-4.84E-9*(PTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11
   ! 
-  ZO10(J) = 0.011*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
+  ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
   ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
   ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2
   ZCT10(J)= ZCH10/SQRT(ZCD10(J))
@@ -280,8 +286,6 @@ ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(:)/ZL10(:)))
 ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(:)/ZL10(:)))
 ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(:)/ZL10(:)))
 !
-ZCHARN(:) = 0.011  
-!
 ZZL(:) = 0.0
 !
 DO J=1,SIZE(PTA)
@@ -293,7 +297,7 @@ DO J=1,SIZE(PTA)
   ENDIF
   !
   !then modify Charnork for high wind speeds Chris Fairall's data
-  IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18-10)
+  IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18.-10.)
   IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
   !
   !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 
@@ -359,7 +363,7 @@ DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
         ZWG(J) = 0.2
       ENDIF
     ENDIF  
-    ZDUWG(J) = SQRT(ZVMOD(J)**2. + ZWG(J)**2.)
+    ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2)
     !
   ENDDO
   !
@@ -379,11 +383,8 @@ ZRF(:)   = 0.
 !
 DO J=1,SIZE(PTA)
   !
-  !            4. 0  roughness length over the ocean PZOSEA
   !
-  IF(ZUSR(J)/=0.0) PZ0SEA(J) = 10./EXP(XKARMAN*ZU10(J)/ZUSR(J))
-  !
-  !            4. 1 transfert coefficients PCD, PCH and PCE 
+  !            4. transfert coefficients PCD, PCH and PCE 
   !                 and neutral PCDN, ZCHN, ZCEN 
   !
   PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
@@ -473,7 +474,12 @@ ZDIRCOSZW(:) = 1.
 ZAC(:) = PCH(:)*ZVMOD(:)
 PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON)
 !
+!       5.3 Z0 and Z0H over sea
+!
+PZ0SEA(:) =  ZCHARN(:) * ZUSTAR2(:) / XG + XVZ0CM * PCD(:) / PCDN(:)
+!
 PZ0HSEA(:) = PZ0SEA(:)
+!
 IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE)
 !
 !-------------------------------------------------------------------------------
diff --git a/src/SURFEX/compute_isba_parameters.F90 b/src/SURFEX/compute_isba_parameters.F90
index e33791f5b..6b78b8a20 100644
--- a/src/SURFEX/compute_isba_parameters.F90
+++ b/src/SURFEX/compute_isba_parameters.F90
@@ -278,12 +278,12 @@ IDECADE2 = IDECADE
                         TPREAP=TREAP,PWATSUP=XWATSUP,PIRRIG=XIRRIG             )
 !
 IF(CISBA=='DIF')THEN
-  IDIM_FULL = SIZE(NWG_LAYER_TOT,1)
+  IF (NPROC>1 .OR. NBLOCKTOT>1) THEN  
+    IDIM_FULL = SIZE(NWG_LAYER_TOT,1)
 !$OMP SINGLE
-  DEALLOCATE(NWG_LAYER_TOT)
-  ALLOCATE(NWG_LAYER_TOT(IDIM_FULL,SIZE(NWG_LAYER,2)))
-!$OMP END SINGLE    
-  IF (NPROC>1 .OR. NBLOCKTOT>1) THEN      
+    DEALLOCATE(NWG_LAYER_TOT)
+    ALLOCATE(NWG_LAYER_TOT(IDIM_FULL,SIZE(NWG_LAYER,2)))
+!$OMP END SINGLE             
     DO JL = 1,SIZE(NWG_LAYER,2)
       IF (HPROGRAM=='ASCII ') THEN
 #ifdef ASC
@@ -303,6 +303,8 @@ IF(CISBA=='DIF')THEN
       ENDIF
     ENDDO
   ELSE
+    DEALLOCATE(NWG_LAYER_TOT)
+    ALLOCATE(NWG_LAYER_TOT(SIZE(NWG_LAYER,1),SIZE(NWG_LAYER,2)))          
     NWG_LAYER_TOT = NWG_LAYER
   ENDIF
   NWG_SIZE = 0
diff --git a/src/SURFEX/ecoclimap2_lai.F90 b/src/SURFEX/ecoclimap2_lai.F90
index a537fcc7e..e6dec3d47 100644
--- a/src/SURFEX/ecoclimap2_lai.F90
+++ b/src/SURFEX/ecoclimap2_lai.F90
@@ -42,7 +42,6 @@ USE MODD_SURF_PAR,       ONLY : XUNDEF
 !
 USE MODD_DATA_COVER_n,   ONLY : XDATA_VEGTYPE, NYEAR
 USE MODD_DATA_COVER,     ONLY : XDATA_LAI, XDATA_LAI_ALL_YEARS, LCLIM_LAI, &
-
                                   NECO2_START_YEAR, NECO2_END_YEAR  
 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE, JPCOVER
 
diff --git a/src/SURFEX/init_isban.F90 b/src/SURFEX/init_isban.F90
index 852791c64..2bb3f96a0 100644
--- a/src/SURFEX/init_isban.F90
+++ b/src/SURFEX/init_isban.F90
@@ -52,6 +52,7 @@ SUBROUTINE INIT_ISBA_n    (HPROGRAM,HINIT,OLAND_USE,                    &
 !!      A.L. Gibelin   06/09 : soil carbon initialisation
 !!      B. Decharme    07/11 : read pgd+prep
 !!      R. Alkama      05/12 : new carbon spinup
+!!      J.Escobar      11/13 : add USE MODI_DEFAULT_CROCUS
 !!
 !-------------------------------------------------------------------------------
 !
diff --git a/src/SURFEX/isba_snow_agr.F90 b/src/SURFEX/isba_snow_agr.F90
index f6cca219b..0e3aba595 100644
--- a/src/SURFEX/isba_snow_agr.F90
+++ b/src/SURFEX/isba_snow_agr.F90
@@ -185,7 +185,7 @@ REAL, DIMENSION(:), INTENT(OUT) :: PGFLUX_ISBA! flux through the ground
 REAL, DIMENSION(:), INTENT(IN)    :: PFFG,PFFV,PFF
 REAL, DIMENSION(:), INTENT(INOUT) :: PLE_FLOOD, PLEI_FLOOD ! Flood evaporation
 !
-REAL, DIMENSION(:),   INTENT(OUT) :: PRI       ! Total Ridcharson number
+REAL, DIMENSION(:),   INTENT(INOUT) :: PRI       ! Total Ridcharson number
 !
 !*      0.2    declarations of local variables
 !
diff --git a/src/SURFEX/modd_surfex_omp.F90 b/src/SURFEX/modd_surfex_omp.F90
index d6af43e2f..bdc955898 100644
--- a/src/SURFEX/modd_surfex_omp.F90
+++ b/src/SURFEX/modd_surfex_omp.F90
@@ -26,6 +26,8 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original       26/06/12
+!!      Modified    11/2013 by J.Escobar :add !$ to inhibit completly omp dependency  
+!!-------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
 !             ------------
@@ -40,7 +42,7 @@ USE PARKIND1  ,ONLY : JPRB
 IMPLICIT NONE
 !
 #ifndef AIX64
-!$INCLUDE 'omp_lib.h'
+!$ INCLUDE 'omp_lib.h'
 #endif
 !
 INTEGER :: NBLOCKTOT = 1
diff --git a/src/SURFEX/mode_psychro.F90 b/src/SURFEX/mode_psychro.F90
index b97b9c360..9c4849116 100644
--- a/src/SURFEX/mode_psychro.F90
+++ b/src/SURFEX/mode_psychro.F90
@@ -27,7 +27,8 @@ MODULE MODE_PSYCHRO
 !!
 !!    MODIFICATIONS
 !!    -------------
-!!      Original    12/04/11 
+!!      Original    12/04/11
+!!      J.Escobar   11/13 :  remove space in ELSEWHERE statement
 !
 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
 USE PARKIND1  ,ONLY : JPRB
diff --git a/src/SURFEX/mode_read_surf_asc.F90 b/src/SURFEX/mode_read_surf_asc.F90
index 807b4d2c5..f95aa49c4 100644
--- a/src/SURFEX/mode_read_surf_asc.F90
+++ b/src/SURFEX/mode_read_surf_asc.F90
@@ -550,7 +550,7 @@ IF (NRANK==NPIO) THEN
   !
 !$OMP END SINGLE
   !
-ELSE
+ELSEIF (HDIR/='-') THEN
 !$OMP SINGLE
   ALLOCATE(NWORKD(0))
 !$OMP END SINGLE
diff --git a/src/SURFEX/mode_read_surf_fa.F90 b/src/SURFEX/mode_read_surf_fa.F90
index f3c8cb3c6..38de166f9 100644
--- a/src/SURFEX/mode_read_surf_fa.F90
+++ b/src/SURFEX/mode_read_surf_fa.F90
@@ -200,11 +200,7 @@ IF (NRANK==NPIO) THEN
   !
 !$OMP END SINGLE
   !
-ELSEIF (HDIR=='-') THEN
-!$OMP SINGLE
-  ALLOCATE(XWORKD(KL))
-!$OMP END SINGLE
-ELSE
+ELSEIF (HDIR/='-') THEN
 !$OMP SINGLE
   ALLOCATE(XWORKD(0))
 !$OMP END SINGLE
@@ -228,10 +224,12 @@ IF (HDIR=='A') THEN  ! no distribution on other tasks
 #endif    
   ENDIF
 ELSEIF (HDIR=='-') THEN ! distribution of the total field on other tasks
-#ifndef NOMPI          
+#ifndef NOMPI
   IF (NPROC>1) THEN
 !$OMP SINGLE
     XTIME0 = MPI_WTIME()
+    CALL MPI_BCAST(NFULL,KIND(NFULL)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
+    IF ( NRANK/=NPIO ) ALLOCATE(XWORKD(NFULL))
     CALL MPI_BCAST(XWORKD(1:KL),KL*KIND(XWORKD)/4,MPI_REAL,NPIO,NCOMM,INFOMPI)
     XTIME_COMM_READ = XTIME_COMM_READ + (MPI_WTIME() - XTIME0)
 !$OMP END SINGLE
@@ -320,7 +318,7 @@ IF (NRANK==NPIO) THEN
   !
 !$OMP SINGLE
   !
-  ALLOCATE(XWORKD2(NFULL,KL2))
+  ALLOCATE(XWORKD2(NFULL,KL2)) 
   !  
   DO JL=1,KL2
     WRITE(YPATCH,'(I2.2)')JL
@@ -336,11 +334,7 @@ IF (NRANK==NPIO) THEN
   !
 !$OMP END SINGLE
   !
-ELSEIF (HDIR=='-') THEN
-!$OMP SINGLE
-  ALLOCATE(XWORKD2(KL1,KL2))
-!$OMP END SINGLE
-ELSE
+ELSEIF (HDIR/='-') THEN
 !$OMP SINGLE
   ALLOCATE(XWORKD2(0,0))
 !$OMP END SINGLE
@@ -368,6 +362,8 @@ ELSEIF (HDIR=='-') THEN ! distribution of the total field on other tasks
 #ifndef NOMPI
   IF (NPROC>1) THEN
     XTIME0 = MPI_WTIME()
+    CALL MPI_BCAST(NFULL,KIND(NFULL)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
+    IF ( NRANK/=NPIO ) ALLOCATE(XWORKD2(NFULL,KL2))        
     CALL MPI_BCAST(XWORKD2(1:KL1,1:KL2),KL1*KL2*KIND(XWORKD2)/4,MPI_REAL,NPIO,NCOMM,INFOMPI)
     XTIME_COMM_READ = XTIME_COMM_READ + (MPI_WTIME() - XTIME0)
   ENDIF
@@ -502,38 +498,36 @@ NWORKB = 0
 XTIME0 = MPI_WTIME()
 #endif
 !
+IF (HDIR=='-') THEN
+!$OMP SINGLE
+  ALLOCATE(NWORKD(KL))
+!$OMP END SINGLE
+ENDIF
+!
 IF (NRANK==NPIO) THEN
   !
 !$OMP SINGLE
-  !
+  !  
   YMASK = CMASK
   CALL IO_BUFF_n(HREC,'R',LWORK0)
   IF (LWORK0) YMASK = 'FULL  '
   !
-  IF (HDIR/='-') THEN
-     ALLOCATE(NWORKD(NFULL))
-  ELSE
-     ALLOCATE(NWORKD(KL))
-  ENDIF
+  IF (HDIR=='A') THEN
+    ALLOCATE(NWORKD(KL))
+  ELSEIF (HDIR/='-') THEN
+    ALLOCATE(NWORKD(NFULL))
+  END IF    
   !
   YNAME = TRIM(YMASK)//TRIM(HREC)
-  IF (HDIR=="-") THEN
-    CALL FALIT_I_D(NWORKB,NUNIT_FA,YNAME,KL,NWORKD(1:KL))
-    IF (NWORKB/=0) CALL ERROR_READ_SURF_FA(HREC,NWORKB)
-  ELSE
-    CALL FALIT_I_D(NWORKB,NUNIT_FA,YNAME,NFULL,NWORKD)
-    IF (NWORKB/=0) CALL ERROR_READ_SURF_FA(HREC,NWORKB)
-  ENDIF
+  !
+  CALL FALIT_I_D(NWORKB,NUNIT_FA,YNAME,SIZE(NWORKD),NWORKD)
+  IF (NWORKB/=0) CALL ERROR_READ_SURF_FA(HREC,NWORKB)
   !
   CWORK0 = YNAME
   !
 !$OMP END SINGLE
   !
-ELSEIF (HDIR=='-') THEN
-!$OMP SINGLE
-  ALLOCATE(NWORKD(KL))
-!$OMP END SINGLE
-ELSE
+ELSEIF (HDIR/='-') THEN
 !$OMP SINGLE
   ALLOCATE(NWORKD(0))
 !$OMP END SINGLE
@@ -561,6 +555,8 @@ ELSEIF (HDIR=='-') THEN ! distribution of the total field on other tasks
 #ifndef NOMPI
   IF (NPROC>1) THEN
     XTIME0 = MPI_WTIME()
+    CALL MPI_BCAST(NFULL,KIND(NFULL)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
+    IF ( NRANK/=NPIO ) ALLOCATE(NWORKD(NFULL))    
     CALL MPI_BCAST(NWORKD(1:KL),KL*KIND(NWORKD)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
     XTIME_COMM_READ = XTIME_COMM_READ + (MPI_WTIME() - XTIME0)
   ENDIF
diff --git a/src/SURFEX/mode_soil.F90 b/src/SURFEX/mode_soil.F90
index 45ca2a0c1..324107ef7 100644
--- a/src/SURFEX/mode_soil.F90
+++ b/src/SURFEX/mode_soil.F90
@@ -36,6 +36,8 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original        3/12/98
+!!
+!!       02/2014       B. Decharme correction for W33 functions
 !-----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -196,8 +198,8 @@ CONTAINS
    SELECT CASE (HPEDOTF)
      CASE ('CH78')
        WHERE(PCLAY/=XUNDEF) &
-         W33_FUNC_2D = 0.3005674207-0.4725429691*PSAND+0.2080416963*PCLAY    &
-                     + 0.3055907979*PSAND**(1./3.)-0.1202046909*PCLAY**(1./3.)
+         W33_FUNC_2D = 0.2298915119-0.4062575773*PSAND+0.0874218705*PCLAY    &
+                     + 0.2942558675*PSAND**(1./3.)+0.0413771051*PCLAY**(1./3.)
      CASE ('CO84')
        WHERE(PCLAY/=XUNDEF) &
          W33_FUNC_2D = 0.2016592588-0.5785747196*PSAND+0.1113006987*PCLAY    &
@@ -544,8 +546,8 @@ CONTAINS
    SELECT CASE (HPEDOTF)
      CASE ('CH78')
        WHERE(PCLAY/=XUNDEF) &
-         W33_FUNC_1D = 0.3005674207-0.4725429691*PSAND+0.2080416963*PCLAY    &
-                     + 0.3055907979*PSAND**(1./3.)-0.1202046909*PCLAY**(1./3.)
+         W33_FUNC_1D = 0.2298915119-0.4062575773*PSAND+0.0874218705*PCLAY    &
+                     + 0.2942558675*PSAND**(1./3.)+0.0413771051*PCLAY**(1./3.)                     
      CASE ('CO84')
        WHERE(PCLAY/=XUNDEF) &
          W33_FUNC_1D = 0.2016592588-0.5785747196*PSAND+0.1113006987*PCLAY    &
diff --git a/src/SURFEX/pgd_surf_atm.F90 b/src/SURFEX/pgd_surf_atm.F90
index a5c57f127..6a5ca9499 100644
--- a/src/SURFEX/pgd_surf_atm.F90
+++ b/src/SURFEX/pgd_surf_atm.F90
@@ -34,6 +34,7 @@
 !!
 !!    Original     13/10/03
 !!      A. Lemonsu      05/2009         Ajout de la clef LGARDEN pour TEB
+!!      J. Escobar      11/2013         Add USE MODI_READ_NAM_PGD_CHEMISTRY
 !----------------------------------------------------------------------------
 !
 !*    0.     DECLARATION
diff --git a/src/SURFEX/pgd_teb_veg.F90 b/src/SURFEX/pgd_teb_veg.F90
index 0ec7cc922..921eaebaa 100644
--- a/src/SURFEX/pgd_teb_veg.F90
+++ b/src/SURFEX/pgd_teb_veg.F90
@@ -34,6 +34,7 @@
 !!    ------------
 !!
 !!    Original    03/2010
+!!    J.Escobar   11/2013   Add USE MODI_PGD_TEB_GREENROOF
 !!
 !----------------------------------------------------------------------------
 !
diff --git a/src/SURFEX/prep_isba_netcdf.F90 b/src/SURFEX/prep_isba_netcdf.F90
index db6365ef2..b69653145 100644
--- a/src/SURFEX/prep_isba_netcdf.F90
+++ b/src/SURFEX/prep_isba_netcdf.F90
@@ -25,6 +25,7 @@ SUBROUTINE PREP_ISBA_NETCDF(HPROGRAM,HSURF,HFILE,KLUOUT,PFIELD)
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    04/2012
+!!      J.Escobar   11/2013  Add USE MODI_GET_TYPE_DIM_n
 !!------------------------------------------------------------------
 !
 
@@ -38,7 +39,7 @@ USE MODE_READ_CDF
 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
 USE PARKIND1  ,ONLY : JPRB
 !
-USE MODI_GET_TYPE_DIM_N
+USE MODI_GET_TYPE_DIM_n
 !
 IMPLICIT NONE
 
diff --git a/src/SURFEX/prep_ocean_ascllv.F90 b/src/SURFEX/prep_ocean_ascllv.F90
index 24719128a..fb8be4fbf 100644
--- a/src/SURFEX/prep_ocean_ascllv.F90
+++ b/src/SURFEX/prep_ocean_ascllv.F90
@@ -36,6 +36,7 @@ SUBROUTINE PREP_OCEAN_ASCLLV(HPROGRAM,HSURF,HFILE, &
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    01/2011
+!!      J.Escobar   11/2013   Add USE MODI_ABOR1_SFX and USE MODI_GET_SURF_MASK_N
 !!------------------------------------------------------------------
 !
 USE MODD_PREP,       ONLY : CINTERP_TYPE, CINGRID_TYPE
@@ -56,7 +57,7 @@ USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
 USE PARKIND1  ,ONLY : JPRB
 !
 USE MODI_ABOR1_SFX
-USE MODI_GET_SURF_MASK_N
+USE MODI_GET_SURF_MASK_n
 !
 IMPLICIT NONE
 !
diff --git a/src/SURFEX/prep_snow_extern.F90 b/src/SURFEX/prep_snow_extern.F90
index 841133f0e..e3254439b 100644
--- a/src/SURFEX/prep_snow_extern.F90
+++ b/src/SURFEX/prep_snow_extern.F90
@@ -5,7 +5,39 @@
 !     #########
 SUBROUTINE PREP_SNOW_EXTERN(HPROGRAM,HSURF,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
                             KLUOUT,PFIELD,OSNOW_IDEAL,KLAYER)
-!     #################################################################################
+
+!!     ##########################################################################
+!
+!
+!!****  *PREP_SNOW_EXTERN*  
+!!
+!!    PURPOSE
+!!    -------
+!       Read and prepare initial snow fields from external files
+!     
+!!**  METHOD
+!!    ------
+!
+!!    EXTERNAL
+!!    --------
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------ 
+!!
+!!      
+!!    REFERENCE
+!!    ---------
+!!
+!!    
+!!    AUTHOR
+!!    ------
+!!	   * Meteo-France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    ?
+!!       02/2014 E. Martin : cor. for passing from from multilayer to a single layer
+!-------------------------------------------------------------------------------
 !
 !
 USE MODD_TYPE_SNOW
@@ -161,10 +193,8 @@ ELSE
   CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
 ENDIF
 !
-IF (TZSNOW%NLAYER.GT.KLAYER) THEN
-  TZSNOW%NLAYER=KLAYER
-ELSEIF (TZSNOW%NLAYER.LT.KLAYER) THEN
-  CALL ABOR1_SFX("PREP_SNOW_EXTERN: SNOW NLAYER IN EXTERN FILE MUST BE GROWER THAN CURRENT NLAYER")
+IF (TZSNOW%NLAYER.LT.KLAYER) THEN
+  CALL ABOR1_SFX("PREP_SNOW_EXTERN: SNOW NLAYER IN EXTERN FILE MUST BE HIGHER THAN CURRENT NLAYER")
 ENDIF
 !
 !
@@ -185,7 +215,8 @@ SELECT CASE (HSURF(1:3))
       ZFIELD = 0.
       DO JLAYER=1,TZSNOW%NLAYER
         ZFIELD(:,1,:) = ZFIELD(:,1,:) + TZSNOW%WSNOW(:,JLAYER,:)
-      END DO
+      END DO 
+      WHERE ( ZFIELD(:,1,:)>XUNDEF ) ZFIELD(:,1,:)=XUNDEF
       ALLOCATE(PFIELD(INI,1,IVEGTYPE))
       CALL PUT_ON_ALL_VEGTYPES(INI,1,IPATCH,IVEGTYPE,ZFIELD,PFIELD)
     ENDIF
diff --git a/src/SURFEX/prep_ver_snow.F90 b/src/SURFEX/prep_ver_snow.F90
index 469b57f9a..0b4601b7a 100644
--- a/src/SURFEX/prep_ver_snow.F90
+++ b/src/SURFEX/prep_ver_snow.F90
@@ -26,6 +26,7 @@ SUBROUTINE PREP_VER_SNOW(TPSNOW,PZS_LS,PZS,PTG_LS,PTG,KDEEP_SOIL)
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    01/2004
+!!       02/2014    E. Martin vertical correction applied to snow cover and not by layers
 !!------------------------------------------------------------------
 !
 
@@ -128,7 +129,7 @@ IF (PRESENT(PTG)) THEN
   DO JPATCH=1,IPATCH
     DO JLAYER=1,TPSNOW%NLAYER
       WHERE(ZWSNOW_LS(:,JLAYER,JPATCH)>0..AND.((PTG(:,KDEEP_SOIL,JPATCH)-XTT >= 2.).OR.(PZS(:) > PZS_LS(:))))
-        ZWSNOW(:,JLAYER,JPATCH) = ZWSNOW_LS(:,JLAYER,JPATCH) + XWSNOW_CLIM_GRAD  * (PZS(:) - PZS_LS(:))  
+        ZWSNOW(:,JLAYER,JPATCH) = ZWSNOW_LS(:,JLAYER,JPATCH) + ( XWSNOW_CLIM_GRAD  * (PZS(:) - PZS_LS(:))/TPSNOW%NLAYER)  
         ZWSNOW(:,JLAYER,JPATCH) = MAX(ZWSNOW(:,JLAYER,JPATCH),0.)
       END WHERE
     END DO
@@ -137,7 +138,7 @@ ELSE
   DO JPATCH=1,IPATCH
     DO JLAYER=1,TPSNOW%NLAYER
       WHERE(ZWSNOW_LS(:,JLAYER,JPATCH)>0.)
-        ZWSNOW(:,JLAYER,JPATCH) = ZWSNOW_LS(:,JLAYER,JPATCH) + XWSNOW_CLIM_GRAD  * (PZS(:) - PZS_LS(:))
+        ZWSNOW(:,JLAYER,JPATCH) = ZWSNOW_LS(:,JLAYER,JPATCH) + ( XWSNOW_CLIM_GRAD  * (PZS(:) - PZS_LS(:))/TPSNOW%NLAYER)
         ZWSNOW(:,JLAYER,JPATCH) = MAX(ZWSNOW(:,JLAYER,JPATCH),0.)
       END WHERE
     END DO
@@ -177,7 +178,7 @@ IF (PRESENT(PTG)) THEN
   ALLOCATE(ZTSNOW2(SIZE(TPSNOW%WSNOW,1),TPSNOW%NLAYER,IPATCH))
   DO JPATCH=1,IPATCH
     DO JLAYER=1,TPSNOW%NLAYER
-      ZWSNOW2(:,JLAYER,JPATCH) =  XWSNOW_CLIM_GRAD  * (PZS(:) - ZZSFREEZE(:,JPATCH))
+      ZWSNOW2(:,JLAYER,JPATCH) =  XWSNOW_CLIM_GRAD  * (PZS(:) - ZZSFREEZE(:,JPATCH))/TPSNOW%NLAYER
       ZWSNOW2(:,JLAYER,JPATCH) = MAX(ZWSNOW2(:,JLAYER,JPATCH),0.)
       ZTSNOW2      (:,JLAYER,JPATCH) = PTG(:,KDEEP_SOIL,JPATCH)
     END DO
diff --git a/src/SURFEX/prep_ver_teb.F90 b/src/SURFEX/prep_ver_teb.F90
index 02a4f0638..527a3346b 100644
--- a/src/SURFEX/prep_ver_teb.F90
+++ b/src/SURFEX/prep_ver_teb.F90
@@ -209,12 +209,12 @@ XT_CANYON = XT_CANYON  + XT_CLIM_GRAD  * (XZS - XZS_LS)
 !* estimation of pressure at large-scale orography
 !
 ALLOCATE(ZP_LS(SIZE(XQ_CANYON)))
-ZP_LS = XP00 * EXP(-(XG/XRD/ZT0)*XZS_LS+(XG*XT_CLIM_GRAD/(2.*XRD*ZT0))*XZS_LS**2)
+ZP_LS = XP00 * EXP(-(XG/XRD/ZT0)*XZS_LS +(XG*XT_CLIM_GRAD/(2.*XRD*ZT0**2))*XZS_LS**2)
 !
 !* estimation of pressure at output orography
 !
 ALLOCATE(ZP(SIZE(XQ_CANYON)))
-ZP    = XP00 * EXP(-(XG/XRD/ZT0)*XZS   +(XG*XT_CLIM_GRAD/(2.*XRD*ZT0))*XZS   **2)
+ZP    = XP00 * EXP(-(XG/XRD/ZT0)*XZS   +(XG*XT_CLIM_GRAD/(2.*XRD*ZT0**2))*XZS   **2)
 !
 !* conservation of estimated relative humidity
 !
diff --git a/src/SURFEX/read_gridtype_gauss.F90 b/src/SURFEX/read_gridtype_gauss.F90
index 8b538302d..56cea6c31 100644
--- a/src/SURFEX/read_gridtype_gauss.F90
+++ b/src/SURFEX/read_gridtype_gauss.F90
@@ -101,7 +101,11 @@ IF (LHOOK) CALL DR_HOOK('READ_GRIDTYPE_GAUSS',0,ZHOOK_HANDLE)
 !
  CALL READ_SURF(HPROGRAM,'NLATI',INLATI,KRESP,HDIR=HDIR)
 ALLOCATE(INLOPA(INLATI))
- CALL READ_SURF(HPROGRAM,'NLOPA',INLOPA(:),KRESP,HDIR='-')
+IF (HDIR=='A') THEN
+  CALL READ_SURF(HPROGRAM,'NLOPA',INLOPA(:),KRESP,HDIR=HDIR)
+ELSE
+  CALL READ_SURF(HPROGRAM,'NLOPA',INLOPA(:),KRESP,HDIR='-')
+ENDIF
  CALL READ_SURF(HPROGRAM,'LATGAUSS',ZLAT(:),KRESP,HDIR=HDIR)
  CALL READ_SURF(HPROGRAM,'LONGAUSS',ZLON(:),KRESP,HDIR=HDIR)
  CALL READ_SURF(HPROGRAM,'LAT_G_XY',ZLAT_XY(:),KRESP,HDIR=HDIR)
diff --git a/src/SURFEX/read_nam_grid_gauss.F90 b/src/SURFEX/read_nam_grid_gauss.F90
index 58708334a..639ee9647 100644
--- a/src/SURFEX/read_nam_grid_gauss.F90
+++ b/src/SURFEX/read_nam_grid_gauss.F90
@@ -33,6 +33,7 @@ SUBROUTINE READ_NAM_GRID_GAUSS(HPROGRAM,KGRID_PAR,KL,PGRID_PAR)
 !!    -------------
 !!      Original    01/2004
 !!      B. Decharme    2008  Comput and save the Mesh size
+!!                     2013  Bug lat and lon for non rotat-strech grid
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -184,9 +185,8 @@ ENDIF
 !
 ALLOCATE(ZLAT_XY(KL))
 ALLOCATE(ZLON_XY(KL))
-
+!
  CALL COMP_GRIDTYPE_GAUSS(INLATI,INLOPA,KL,ITYP,ZLAT_XY,ZLON_XY)
-
 !
 !---------------------------------------------------------------------------
 !
@@ -196,7 +196,13 @@ ALLOCATE(ZLON_XY(KL))
 !* all points are used
 ALLOCATE(ZLAT(KL))
 ALLOCATE(ZLON(KL))
- CALL LATLON_GAUSS(ZLON_XY,ZLAT_XY,KL,ZLOPO,ZLAPO,ZCODIL,ZLON,ZLAT)
+!
+IF(ZCODIL==1.0.AND.ZLAPO==90.0.AND.ZLOPO==0.0)THEN
+  ZLON(:)=ZLON_XY(:)
+  ZLAT(:)=ZLAT_XY(:)
+ELSE
+  CALL LATLON_GAUSS(ZLON_XY,ZLAT_XY,KL,ZLOPO,ZLAPO,ZCODIL,ZLON,ZLAT)
+ENDIF
 !
 !---------------------------------------------------------------------------
 !
diff --git a/src/SURFEX/snowcro.F90 b/src/SURFEX/snowcro.F90
index f2edd50ea..b823a34c5 100644
--- a/src/SURFEX/snowcro.F90
+++ b/src/SURFEX/snowcro.F90
@@ -3450,7 +3450,7 @@ DO JJ=1, SIZE(PSNOW(:))
       !remaining snow for remaining layers
       THICKNESS_INTERMEDIATE=PSNOW_UPPER - SUM(ZDZOPT(JJ,1:5))-ZDZOPT(JJ,NB_UPPER_LAYER)
 
-      IF (PSNOW_UPPER<DEPTH_THRESHOLD1) THEN
+      IF ((PSNOW_UPPER<=DEPTH_THRESHOLD1) .OR.(NB_UPPER_LAYER<8)) THEN
         NB_INTERMEDIATE=NB_UPPER_LAYER-6
         END_INTERMEDIATE=NB_UPPER_LAYER-1
       ELSE
@@ -3940,42 +3940,36 @@ DO JJ=1,SIZE(PSNOW(:))
       INLVLS_USE(JJ)=INLVLS_USE(JJ)-1
       PSNOWDZN(JJ,INLVLS_USE(JJ))=PSNOWDZN(JJ,INLVLS_USE(JJ))+  &
                    PSNOWDZN(JJ,INLVLS_USE(JJ)+1)
+      PSNOWDZN(JJ,INLVLS_USE(JJ)+1) = 0.
    ELSE
 ! internal layer
-       DO JST=INLVLS_USE(JJ)-1,4 
+       DO JST=INLVLS_USE(JJ)-1,4,-1 
           IF(ABS(PSNOWDZN(JJ,JST)-PSNOWDZ(JJ,JST)) &
               > UEPSI) EXIT ! old/new grid differ ==> go to next grid point 
-         ZDIFTYPE_SUP= SNOW3LDIFTYP(PSNOWGRAN1(JJ,JST-1),PSNOWGRAN1(JJ, JST),&
-                       PSNOWGRAN2(JJ,JST-1),PSNOWGRAN2(JJ, JST))
-         ZDIFTYPE_SUP=MAX(DIFF_1,MIN(DIFF_MAX,ZDIFTYPE_SUP))
-         IF(PSNOWDZN(JJ,JST)>ZDZOPT(JJ,JST)*AGREG_COEF_1/ZDIFTYPE_SUP &
-            .OR.PSNOWDZN(JJ,JST)+PSNOWDZN(JJ,JST-1)>AGREG_COEF_2*&
-             MAX(ZDZOPT(JJ,JST),ZDZOPT(JJ,JST-1))) CYCLE! cheks upper layer
-! shallow internal layer to be merged with the upper layer             
-         PSNOWDZN(JJ,JST)= PSNOWDZN(JJ,JST)+PSNOWDZN(JJ,JST-1)
+         IF(PSNOWDZN(JJ,JST)> 0.001) CYCLE! 
+! If an internal layer is too shallow, it is merged with the upper layer             
+         PSNOWDZN(JJ,JST-1)= PSNOWDZN(JJ,JST)+PSNOWDZN(JJ,JST-1)
          INLVLS_USE(JJ)=INLVLS_USE(JJ)-1
-! shifts the upper layers       
-         DO JST_1=JST-1,1
-            PSNOWDZN(JJ,JST_1)=PSNOWDZ(JJ,JST_1-1)
-            ZDZOPT(JJ,JST_1)=ZDZOPT(JJ,JST_1-1)
-         ENDDO             
+! shifts the lower layers       
+         DO JST_1=JST,INLVLS_USE(JJ)
+            PSNOWDZN(JJ,JST_1)=PSNOWDZ(JJ,JST_1+1)
+            ZDZOPT(JJ,JST_1)=ZDZOPT(JJ,JST_1+1)
+         ENDDO
+         PSNOWDZN(JJ,INLVLS_USE(JJ)+1) = 0.
          EXIT ! goto to next grid point
       ENDDO ! end loop internal layers             
    ENDIF
 ENDDO ! end grid loops for checking shallow layers   
-
-
-
-
 !
 !final check of the consistensy of the new grid size
 !
 DO JJ=1,SIZE(PSNOW(:))
   IF(ABS(SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ)))-PSNOW(JJ)) > UEPSI) THEN
-   write(*,*) 'error in grid resizing',JJ, SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ))), &
-               PSNOW(JJ),INLVLS_USE(JJ) 
-write( *,*) 'PSNOWDZ:',     PSNOWDZ
-write( *,*) 'PSNOWDZN:',     PSNOWDZN
+   write(*,*) 'error in grid resizing JJ,USE,H_new,H_old,delta,ZSNOWFALL',&
+            JJ,INLVLS_USE(JJ), SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ))), PSNOW(JJ), & 
+        SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ)))-PSNOW(JJ),ZSNOWFALL(JJ)
+   write( *,*) 'JJ , PSNOWDZ(JJ):',JJ ,  PSNOWDZ(JJ,:)
+   write( *,*) 'JJ , PSNOWDZN(JJ):',JJ , PSNOWDZN(JJ,:)
  CALL ABOR1_SFX("SNOWCRO: error in grid resizing")
    ENDIF       
 ENDDO
diff --git a/src/SURFEX/sunpos.F90 b/src/SURFEX/sunpos.F90
index f5ea2be40..e2efdbcdf 100644
--- a/src/SURFEX/sunpos.F90
+++ b/src/SURFEX/sunpos.F90
@@ -48,6 +48,7 @@
 !!      (J.Stein)            01:04/96  bug correction for ZZEANG     
 !!      (K. Suhre)           14/02/97  bug correction for ZLON0     
 !!      (V. Masson)          01/03/03  add zenithal angle output
+!!      (J.Escobar)           11/2013  add !$ to inhibit completly omp dependency  
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -66,7 +67,7 @@ USE PARKIND1  ,ONLY : JPRB
 IMPLICIT NONE
 !
 #ifndef AIX64
-!$INCLUDE 'omp_lib.h'
+!$ INCLUDE 'omp_lib.h'
 #endif
 !
 !*       0.1   Declarations of dummy arguments :
-- 
GitLab