Index: LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 2672)
@@ -303,4 +303,20 @@
        CALL getin('rugos',rugos)
 
+!Config  Key  = rugosh
+!Config  Desc = coefficient de frottement
+!Config  Def  = rugos
+!Config  Help = calcul du Cdrag
+       rugosh = rugos
+       CALL getin('rugosh',rugosh)
+
+
+
+!Config  Key  = snowmass
+!Config  Desc = mass de neige de la surface en kg/m2
+!Config  Def  = 0.0000
+!Config  Help = snowmass
+       snowmass = 0.0000
+       CALL getin('snowmass',snowmass)
+
 !Config  Key  = wtsurf et wqsurf
 !Config  Desc = ???
@@ -395,4 +411,5 @@
       write(lunout,*)' zsurf = ', zsurf
       write(lunout,*)' rugos = ', rugos
+      write(lunout,*)' snowmass=', snowmass
       write(lunout,*)' wtsurf = ', wtsurf
       write(lunout,*)' wqsurf = ', wqsurf
@@ -3141,4 +3158,91 @@
 
 !======================================================================
+        SUBROUTINE interp_gabls4_time(day,day1,annee_ref                              &
+     &             ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4    &
+     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                          &
+     &             ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
+        implicit none
+
+!---------------------------------------------------------------------------------------
+! Time interpolation of a 2D field to the timestep corresponding to day
+!
+! day: current julian day 
+! day1: first day of the simulation
+! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
+! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
+!---------------------------------------------------------------------------------------
+
+#include "compar1d.h"
+
+! inputs:
+        integer annee_ref
+        integer nt_gabls4,nlev_gabls4
+        integer year_ini_gabls4
+        real day, day1,day_ini_gabls4,dt_gabls4
+        real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
+        real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
+        real tg_gabls4(nt_gabls4), tg_prof
+! outputs:
+        real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
+        real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
+! local:
+        integer it_gabls41, it_gabls42,k
+        real timeit,time_gabls41,time_gabls42,frac
+
+
+
+! Check that initial day of the simulation consistent with gabls4 period:
+       if (forcing_type.eq.8 ) then
+       print *,'annee_ref=',annee_ref
+       print *,'day1=',day1
+       print *,'day_ini_gabls4=',day_ini_gabls4
+       if (annee_ref.ne.2009) then
+        print*,'Pour gabls4, annee_ref doit etre 2009'
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
+        print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
+        print*,'Changer dayref dans run.def',day1,day_ini_gabls4
+        stop
+       endif
+       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
+        print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
+        print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
+        stop
+       endif
+       endif
+
+      timeit=(day-day_ini_gabls4)*86400
+       print *,'day,day_ini_gabls4=',day,day_ini_gabls4
+       print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
+
+! Determine the closest observation times:
+       it_gabls41=INT(timeit/dt_gabls4)+1
+       it_gabls42=it_gabls41 + 1
+       time_gabls41=(it_gabls41-1)*dt_gabls4
+       time_gabls42=(it_gabls42-1)*dt_gabls4
+
+       if (it_gabls41 .ge. nt_gabls4) then
+        write(*,*) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
+        stop
+       endif
+
+! time interpolation:
+       frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
+       frac=max(frac,0.0)
+
+
+       do k=1,nlev_gabls4
+        ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
+        vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
+        ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
+        hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
+        enddo
+        tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
+        return
+        END
+
+!======================================================================
         SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
      &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
@@ -4095,4 +4199,308 @@
          end subroutine read_dice
 !=====================================================================
+      subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol                    &
+     &     ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
+
+!program reading initial profils and forcings of the Gabls4 case study
+
+
+      implicit none
+
+#include "netcdf.inc"
+
+      integer ntime,nlevel,nsol
+      integer l,k
+      character*80 :: fich_gabls4
+      real*8 time(ntime)
+
+!  ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees 
+! dans un ordre inverse par rapport a la convention LMDZ
+! ==> il faut tout inverser  (MPL 20141024)
+! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
+      real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
+      real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
+      real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
+
+      real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
+      real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
+      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
+
+      real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
+      real*8 tg(ntime)
+      integer nid, ierr
+      integer nbvar3d
+      parameter(nbvar3d=30)
+      integer var3didin(nbvar3d)
+
+      ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'ERROR: Pb opening forcings nc file '
+         write(*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+
+
+       ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'height'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'depth_sn'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Ug'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Vg'
+      endif
+       ierr=NF_INQ_VARID(nid,"pf",var3didin(5)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'pf'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'theta'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'tempe'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'qv'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"u",var3didin(9))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'u'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"v",var3didin(10))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'v'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvt'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvq'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'tsnow'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'snow_density'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Tg'
+      endif
+
+
+!dimensions lecture
+!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+
+! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
+         do k=1,nlevel
+           zz(k)=zz_i(nlevel+1-k)
+           ug(k,:)=ug_i(nlevel+1-k,:)
+           vg(k,:)=vg_i(nlevel+1-k,:)
+           pf(k)=pf_i(nlevel+1-k)
+           print *,'pf=',pf(k)
+           th(k)=th_i(nlevel+1-k)
+           t(k)=t_i(nlevel+1-k)
+           qv(k)=qv_i(nlevel+1-k)
+           u(k)=u_i(nlevel+1-k)
+           v(k)=v_i(nlevel+1-k)
+           hadvt(k,:)=hadvt_i(nlevel+1-k,:)
+           hadvq(k,:)=hadvq_i(nlevel+1-k,:)
+         enddo
+         return 
+ end subroutine read_gabls4
+!=====================================================================
+
 !     Reads CIRC input files      
 
Index: LMDZ5/trunk/libf/phylmd/dyn1d/1D_decl_cases.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/1D_decl_cases.h	(revision 2672)
@@ -33,5 +33,5 @@
 
         real w_mod(llm), t_mod(llm),q_mod(llm)
-        real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm)
+        real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm),ug_mod(llm),vg_mod(llm)
         real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm)
         real th_mod(llm)
@@ -94,4 +94,38 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!Declarations specifiques au cas GABLS4   (MPL 20141023)
+        character*80 :: fich_gabls4
+        integer nlev_gabls4, nt_gabls4, nsol_gabls4
+        parameter (nlev_gabls4=90, nt_gabls4=37, nsol_gabls4=19)  
+        integer year_ini_gabls4, day_ini_gabls4, mth_ini_gabls4
+        real heure_ini_gabls4
+        real day_ju_ini_gabls4   ! Julian day of gabls4 first day
+        parameter (year_ini_gabls4=2009) 
+        parameter (mth_ini_gabls4=12)
+        parameter (day_ini_gabls4=11)  ! 11 = 11 decembre 2009
+        parameter (heure_ini_gabls4=0.) !0UTC en secondes
+        real dt_gabls4
+        parameter (dt_gabls4=3600.) ! 1 forcage ttes les heures
+
+!profils initiaux:
+        real plev_gabls4(nlev_gabls4)
+        real zz_gabls4(nlev_gabls4)
+        real th_gabls4(nlev_gabls4),t_gabls4(nlev_gabls4),qv_gabls4(nlev_gabls4)
+        real u_gabls4(nlev_gabls4), v_gabls4(nlev_gabls4)
+        real depth_sn_gabls4(nsol_gabls4),tsnow_gabls4(nsol_gabls4),snow_dens_gabls4(nsol_gabls4)
+        real t_gabi(nlev_gabls4),qv_gabi(nlev_gabls4)
+        real u_gabi(nlev_gabls4), v_gabi(nlev_gabls4),ug_gabi(nlev_gabls4), vg_gabi(nlev_gabls4)
+        real ht_gabi(nlev_gabls4),hq_gabi(nlev_gabls4),poub(nlev_gabls4)
+        
+!forcings
+        real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
+        real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
+        real tg_gabls4(nt_gabls4)
+        real ht_profg(nlev_gabls4),hq_profg(nlev_gabls4)
+        real ug_profg(nlev_gabls4),vg_profg(nlev_gabls4)
+        real tg_profg
+         
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 !Declarations specifiques au cas DICE     (MPL 02072013)
         character*80 :: fich_dice
Index: LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 2672)
@@ -192,4 +192,43 @@
       endif ! forcing_dice
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Interpolation gabls4 forcing 
+!---------------------------------------------------------------------
+      if (forcing_gabls4 ) then
+
+       if (prt_level.ge.1) then
+        print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_gabls4=',&
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_gabls4
+       endif
+
+! time interpolation:
+      CALL interp_gabls4_time(daytime,day1,annee_ref                                     &
+     &             ,year_ini_gabls4,day_ju_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4  & 
+     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                            &
+     &             ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg)
+
+        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
+
+! vertical interpolation:
+! on re-utilise le programme interp_dice_vertical: les transformations sur 
+! plev_gabls4,th_gabls4,qv_gabls4,u_gabls4,v_gabls4 ne sont pas prises en compte.
+! seules celles sur ht_profg,hq_profg,ug_profg,vg_profg sont prises en compte.
+
+      CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4         &
+!    &         ,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,poub            &
+     &         ,poub,poub,poub,poub,poub                             &
+     &         ,ht_profg,hq_profg,ug_profg,vg_profg,poub,poub        &
+     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                      &
+     &         ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
+
+      do l = 1, llm
+       ug(l)= ug_mod(l)
+       vg(l)= vg_mod(l)
+       d_th_adv(l)=ht_mod(l)
+       d_q_adv(l,1)=hq_mod(l)
+      enddo
+
+      endif ! forcing_gabls4
+!---------------------------------------------------------------------
+
 !---------------------------------------------------------------------
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 2672)
@@ -473,4 +473,116 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
+! Forcing from GABLS4 experiment
+!---------------------------------------------------------------------
+
+!!!! Si la temperature de surface n'est pas imposée:
+  
+      if (forcing_gabls4) then
+!read GABLS4 forcings
+      
+      fich_gabls4='gabls4_driver.nc'
+      
+       	
+      call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 &
+     & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4)
+
+      write(*,*) 'Forcing GABLS4 lu'
+
+!champs initiaux:
+      do k=1,nlev_gabls4
+         t_gabi(k)=t_gabls4(k)
+         qv_gabi(k)=qv_gabls4(k)
+         u_gabi(k)=u_gabls4(k)
+         v_gabi(k)=v_gabls4(k)
+         poub(k)=0.
+         ht_gabi(k)=ht_gabls4(k,1)
+         hq_gabi(k)=hq_gabls4(k,1)
+         ug_gabi(k)=ug_gabls4(k,1)
+         vg_gabi(k)=vg_gabls4(k,1)
+      enddo 
+  
+      omega(:)=0.      
+      omega2(:)=0.
+      rho(:)=0.
+! vertical interpolation using TOGA interpolation routine:
+!      write(*,*)'avant interp vert', t_proftwp
+!
+!     CALL interp_dice_time(daytime,day1,annee_ref
+!    i             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
+!    i             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
+!    i             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
+!    i             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
+!    o             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
+!    o             ,ustar_prof,psurf_prof,ug_profd,vg_profd
+!    o             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
+!    o             ,omega_profd)
+
+      CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4       &
+     &         ,t_gabi,qv_gabi,u_gabi,v_gabi,poub                  &
+     &         ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub          &
+     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                    &
+     &         ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
+
+! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024
+!     ht_mod(:)=ht_mod(:)/86400.
+!     hq_mod(:)=hq_mod(:)/86400.
+
+! initial and boundary conditions :
+      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
+      do l = 1, llm
+! Ligne du dessous à decommenter si on lit theta au lieu de temp
+!      temp(l) = th_mod(l)*(play(l)/pzero)**rkappa 
+       temp(l) = t_mod(l) 
+       q(l,1) = qv_mod(l)
+       q(l,2) = 0.0
+!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       ug(l)=ug_mod(l)
+       vg(l)=vg_mod(l)
+       
+!
+!	tg=tsurf
+!       
+
+       print *,'***** tsurf=',tsurf
+       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+!      omega(l) = w_mod(l)*(-rg*rho(l))
+       omega(l) = omega_mod(l)
+       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+       
+    
+
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!on applique le forcage total au premier pas de temps
+!attention: signe different de toga
+!      d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
+!forcage en th
+       d_th_adv(l) = ht_mod(l)
+       d_q_adv(l,1) = hq_mod(l)
+       d_q_adv(l,2) = 0.0
+       dt_cooling(l)=0.
+      enddo     
+
+!--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024---------------
+! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par 
+! le coefficient de trainee en surface cd**2=ustar*vent(k=1)
+! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
+! MPL 05082013
+!     ust=ustar_dice(1)
+!     tg=tg_dice(1)
+!     print *,'ust= ',ust
+!     IF (tsurf .LE. 0.) THEN
+!      tsurf= tg_dice(1)
+!     ENDIF
+!     psurf= psurf_dice(1)
+!     solsw_in = (1.-albedo)/albedo*swup_dice(1)
+!     sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1)
+!     PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
+!--------------------------------------------------------------------------------------
+      endif !forcing_gabls4
+
+
+
 ! Forcing from Arm_Cu case                   
 ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
Index: LMDZ5/trunk/libf/phylmd/dyn1d/compar1d.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/compar1d.h	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/compar1d.h	(revision 2672)
@@ -9,4 +9,5 @@
       real :: tsurf
       real :: rugos
+      real :: rugosh
       real :: xqsol(1:2)
       real :: qsurf
@@ -14,4 +15,5 @@
       real :: zsurf
       real :: albedo
+      real :: snowmass
 
       real :: time
@@ -31,10 +33,10 @@
 
       common/com_par1d/                                                 &
-     & nat_surf,tsurf,rugos,                                            &
+     & nat_surf,tsurf,rugos,rugosh,                                     &
      & xqsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,   &
      & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
      & forcing_type,tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo,       &
      & nudge_u,nudge_v,nudge_w,nudge_t,nudge_q,                         &
-     & iflag_nudge,                                                     &
+     & iflag_nudge,snowmass,                                            &
      & restart,ok_old_disvert
 
Index: LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2668)
+++ LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2672)
@@ -21,4 +21,5 @@
        zgam, zmax0, zmea, zpic, zsig, &
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
+  
    USE dimphy
    USE surface_data, only : type_ocean,ok_veget
@@ -131,4 +132,6 @@
         logical :: forcing_amma    = .false.
         logical :: forcing_dice    = .false.
+        logical :: forcing_gabls4  = .false.
+
         logical :: forcing_GCM2SCM = .false.
         logical :: forcing_GCSSold = .false.
@@ -174,5 +177,5 @@
       real :: pzero=1.e5
       real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
-      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
+      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
 
 !---------------------------------------------------------------------
@@ -322,4 +325,6 @@
 !             Different stages: soil model alone, atm. model alone
 !             then both models coupled
+!forcing_type = 8 ==> forcing_gabls4 = .true.
+!             initial profiles and large scale forcings in gabls4_driver.nc
 !forcing_type >= 100 ==> forcing_case = .true.
 !             initial profiles and large scale forcings in cas.nc
@@ -363,4 +368,6 @@
       elseif (forcing_type .eq.7) THEN
        forcing_dice = .true.
+      elseif (forcing_type .eq.8) THEN
+       forcing_gabls4 = .true.
       elseif (forcing_type .eq.101) THEN ! Cindynamo starts 1-10-2011 0h
        forcing_case = .true.
@@ -457,4 +464,5 @@
       IF(forcing_type .EQ. 6) fnday=64800./86400.
 !     IF(forcing_type .EQ. 6) fnday=50400./86400.
+ IF(forcing_type .EQ. 8 ) fnday=129600./86400. 
       annee_ref = anneeref
       mois = 1
@@ -487,4 +495,9 @@
      & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             & 
      & ,day_ju_ini_dice)
+ ELSEIF (forcing_type .eq.8 ) THEN
+! Convert the initial date of GABLS4 to Julian day
+      call ymds2ju                                                         &
+     & (year_ini_gabls4,mth_ini_gabls4,day_ini_gabls4,heure_ini_gabls4     & 
+     & ,day_ju_ini_gabls4)
       ELSEIF (forcing_type .gt.100) THEN
 ! Convert the initial date to Julian day
@@ -699,9 +712,9 @@
 
         fder=0.
-        snsrf(1,:)=0.        ! couverture de neige des sous surface
+        snsrf(1,:)=snowmass ! masse de neige des sous surface
         qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
         fevap=0.
         z0m(1,:)=rugos     ! couverture de neige des sous surface
-        z0h(1,:)=rugos     ! couverture de neige des sous surface
+        z0h(1,:)=rugosh    ! couverture de neige des sous surface
         agesno  = xagesno
         tsoil(:,:,:)=tsurf
@@ -726,11 +739,27 @@
         print*,'avant phyredem'
         pctsrf(1,:)=0.
-        if (nat_surf.eq.0.) then
+          if (nat_surf.eq.0.) then
           pctsrf(1,is_oce)=1.
           pctsrf(1,is_ter)=0.
-        else
+          pctsrf(1,is_lic)=0.
+          pctsrf(1,is_sic)=0.
+        else if (nat_surf .eq. 1) then 
           pctsrf(1,is_oce)=0.
           pctsrf(1,is_ter)=1.
-        end if
+          pctsrf(1,is_lic)=0.
+          pctsrf(1,is_sic)=0.
+        else if (nat_surf .eq. 2) then 
+          pctsrf(1,is_oce)=0.
+          pctsrf(1,is_ter)=0.
+          pctsrf(1,is_lic)=1.
+          pctsrf(1,is_sic)=0.
+        else if (nat_surf .eq. 3) then 
+          pctsrf(1,is_oce)=0.
+          pctsrf(1,is_ter)=0.
+          pctsrf(1,is_lic)=0.
+          pctsrf(1,is_sic)=1.
+
+     end if
+
 
         print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
