Changeset 2720 for LMDZ5/branches/testing/libf/phylmd/dyn1d/1DUTILS.h

Timestamp:

Nov 30, 2016, 1:28:41 PM (8 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2664:2719 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/dyn1d/1DUTILS.h (modified) (21 diffs)

LMDZ5/branches/testing/libf/phylmd/dyn1d/1DUTILS.h

@@ -55 +55 @@
 
 !Config  Key  = prt_level
-!Config  Desc = niveau d'impressions de débogage
+!Config  Desc = niveau d'impressions de debogage
 !Config  Def  = 0
-!Config  Help = Niveau d'impression pour le débogage
+!Config  Help = Niveau d'impression pour le debogage
 !Config         (0 = minimum d'impression)
 !      prt_level = 0
@@ -118 +118 @@
 !             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 
 !             Radiation to be switched off
+!         > 100 ==> forcing_case = .true. or forcing_case2 = .true.
+!             initial profiles from case.nc file 
 !
        forcing_type = 0
@@ -134 +136 @@
         ENDIF
 
-!Paramètres de forçage
+!Parametres de forcage
 !Config  Key  = tend_t
 !Config  Desc = forcage ou non par advection de T
@@ -303 +305 @@
        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 = ???
@@ -342 +360 @@
 !Config  Key  = zpicinp
 !Config  Desc = denivellation orographie
-!Config  Def  = 300.
+!Config  Def  = 0.
 !Config  Help =  input brise
-       zpicinp = 300.
+       zpicinp = 0.
        CALL getin('zpicinp',zpicinp)
 !Config key = nudge_tsoil
@@ -378 +396 @@
        CALL getin('tau_soil_nudge',tau_soil_nudge)
 
+!----------------------------------------------------------
+! Param??tres de for??age pour les forcages communs:
+! Pour les forcages communs: ces entiers valent 0 ou 1
+! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
+! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale
+! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv)
+! forcages en omega, w, vent geostrophique ou ustar
+! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging
+!----------------------------------------------------------
+
+!Config  Key  = tadv
+!Config  Desc = forcage ou non par advection totale de T
+!Config  Def  = false
+!Config  Help = forcage ou non par advection totale de T
+       tadv =0
+       CALL getin('tadv',tadv)
+
+!Config  Key  = tadvv
+!Config  Desc = forcage ou non par advection verticale de T
+!Config  Def  = false
+!Config  Help = forcage ou non par advection verticale de T
+       tadvv =0
+       CALL getin('tadvv',tadvv)
+
+!Config  Key  = tadvh
+!Config  Desc = forcage ou non par advection horizontale de T
+!Config  Def  = false
+!Config  Help = forcage ou non par advection horizontale de T
+       tadvh =0
+       CALL getin('tadvh',tadvh)
+
+!Config  Key  = thadv
+!Config  Desc = forcage ou non par advection totale de Theta
+!Config  Def  = false
+!Config  Help = forcage ou non par advection totale de Theta
+       thadv =0
+       CALL getin('thadv',thadv)
+
+!Config  Key  = thadvv
+!Config  Desc = forcage ou non par advection verticale de Theta
+!Config  Def  = false
+!Config  Help = forcage ou non par advection verticale de Theta
+       thadvv =0
+       CALL getin('thadvv',thadvv)
+
+!Config  Key  = thadvh
+!Config  Desc = forcage ou non par advection horizontale de Theta
+!Config  Def  = false
+!Config  Help = forcage ou non par advection horizontale de Theta
+       thadvh =0
+       CALL getin('thadvh',thadvh)
+
+!Config  Key  = qadv
+!Config  Desc = forcage ou non par advection totale de Q
+!Config  Def  = false
+!Config  Help = forcage ou non par advection totale de Q
+       qadv =0
+       CALL getin('qadv',qadv)
+
+!Config  Key  = qadvv
+!Config  Desc = forcage ou non par advection verticale de Q
+!Config  Def  = false
+!Config  Help = forcage ou non par advection verticale de Q
+       qadvv =0
+       CALL getin('qadvv',qadvv)
+
+!Config  Key  = qadvh
+!Config  Desc = forcage ou non par advection horizontale de Q
+!Config  Def  = false
+!Config  Help = forcage ou non par advection horizontale de Q
+       qadvh =0
+       CALL getin('qadvh',qadvh)
+
+!Config  Key  = trad
+!Config  Desc = forcage ou non par tendance radiative
+!Config  Def  = false
+!Config  Help = forcage ou non par tendance radiative
+       trad =0
+       CALL getin('trad',trad)
+
+!Config  Key  = forc_omega
+!Config  Desc = forcage ou non par omega
+!Config  Def  = false
+!Config  Help = forcage ou non par omega
+       forc_omega =0
+       CALL getin('forc_omega',forc_omega)
+
+!Config  Key  = forc_w
+!Config  Desc = forcage ou non par w
+!Config  Def  = false
+!Config  Help = forcage ou non par w
+       forc_w =0
+       CALL getin('forc_w',forc_w)
+
+!Config  Key  = forc_geo
+!Config  Desc = forcage ou non par geo
+!Config  Def  = false
+!Config  Help = forcage ou non par geo
+       forc_geo =0
+       CALL getin('forc_geo',forc_geo)
+
+! Meme chose que ok_precr_ust
+!Config  Key  = forc_ustar
+!Config  Desc = forcage ou non par ustar
+!Config  Def  = false
+!Config  Help = forcage ou non par ustar
+       forc_ustar =0
+       CALL getin('forc_ustar',forc_ustar)
+       IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
+
+!Config  Key  = nudging_u
+!Config  Desc = forcage ou non par nudging sur u
+!Config  Def  = false
+!Config  Help = forcage ou non par nudging sur u
+       nudging_u =0
+       CALL getin('nudging_u',nudging_u)
+
+!Config  Key  = nudging_v
+!Config  Desc = forcage ou non par nudging sur v
+!Config  Def  = false
+!Config  Help = forcage ou non par nudging sur v
+       nudging_v =0
+       CALL getin('nudging_v',nudging_v)
+
+!Config  Key  = nudging_w
+!Config  Desc = forcage ou non par nudging sur w
+!Config  Def  = false
+!Config  Help = forcage ou non par nudging sur w
+       nudging_w =0
+       CALL getin('nudging_w',nudging_w)
+
+!Config  Key  = nudging_q
+!Config  Desc = forcage ou non par nudging sur q
+!Config  Def  = false
+!Config  Help = forcage ou non par nudging sur q
+       nudging_q =0
+       CALL getin('nudging_q',nudging_q)
+
+!Config  Key  = nudging_t
+!Config  Desc = forcage ou non par nudging sur t
+!Config  Def  = false
+!Config  Help = forcage ou non par nudging sur t
+       nudging_t =0
+       CALL getin('nudging_t',nudging_t)
 
 
@@ -395 +557 @@
       write(lunout,*)' zsurf = ', zsurf
       write(lunout,*)' rugos = ', rugos
+      write(lunout,*)' snowmass=', snowmass
       write(lunout,*)' wtsurf = ', wtsurf
       write(lunout,*)' wqsurf = ', wqsurf
@@ -406 +569 @@
       write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
       write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
+      write(lunout,*)' tadv =      ', tadv
+      write(lunout,*)' tadvv =     ', tadvv
+      write(lunout,*)' tadvh =     ', tadvh
+      write(lunout,*)' thadv =     ', thadv
+      write(lunout,*)' thadvv =    ', thadvv
+      write(lunout,*)' thadvh =    ', thadvh
+      write(lunout,*)' qadv =      ', qadv
+      write(lunout,*)' qadvv =     ', qadvv
+      write(lunout,*)' qadvh =     ', qadvh
+      write(lunout,*)' trad =      ', trad
+      write(lunout,*)' forc_omega = ', forc_omega
+      write(lunout,*)' forc_w     = ', forc_w
+      write(lunout,*)' forc_geo   = ', forc_geo
+      write(lunout,*)' forc_ustar = ', forc_ustar
+      write(lunout,*)' nudging_u  = ', nudging_u
+      write(lunout,*)' nudging_v  = ', nudging_v
+      write(lunout,*)' nudging_t  = ', nudging_t
+      write(lunout,*)' nudging_q  = ', nudging_q
       IF (forcing_type .eq.40) THEN
         write(lunout,*) '--- Forcing type GCSS Old --- with:'
@@ -1106 +1287 @@
 !----------------------------------------------------------------------
 !   Calcul de l'advection verticale (ascendance et subsidence) de 
-!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
-!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou 
+!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
+!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 
 !   sans WTG rajouter une advection horizontale 
 !----------------------------------------------------------------------  
@@ -1180 +1361 @@
 !----------------------------------------------------------------------
 !   Calcul de l'advection verticale (ascendance et subsidence) de 
-!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
-!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou 
+!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
+!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 
 !   sans WTG rajouter une advection horizontale 
 !----------------------------------------------------------------------  
@@ -2934 +3115 @@
        endif
        if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
-        print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
+        print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'
         print*,'Changer dayref dans run.def'
         stop
@@ -3141 +3322 @@
 
 !======================================================================
+        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         &
@@ -3679 +3947 @@
 !=====================================================================
       subroutine read_dice(fich_dice,nlevel,ntime                         &
-     &     ,zz,pres,th,qv,u,v,o3                                          &
+     &     ,zz,pres,t,qv,u,v,o3                                          &
      &     ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg                        &
      &     ,hadvt,hadvq,hadvu,hadvv,w,omega)
@@ -3689 +3957 @@
 
 #include "netcdf.inc"
+#include "YOMCST.h"
 
       integer ntime,nlevel
@@ -3696 +3965 @@
       real*8 zz(nlevel)
 
-      real*8 th(nlevel),pres(nlevel)
+      real*8 th(nlevel),pres(nlevel),t(nlevel)
       real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
       real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
@@ -3702 +3971 @@
       real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
       real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
+      real*8 pzero
 
       integer nid, ierr
@@ -3708 +3978 @@
       integer var3didin(nbvar3d)
 
+      pzero=100000.
       ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
       if (ierr.NE.NF_NOERR) then
@@ -3882 +4153 @@
          endif
 !          write(*,*)'lecture th ok',th
+           do k=1,nlevel
+             t(k)=th(k)*(pres(k)/pzero)**rkappa
+           enddo
 
 #ifdef NC_DOUBLE
@@ -4095 +4369 @@
          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      
 
@@ -4390 +4968 @@
 !
 !  Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
-!   qui n'était pas correcte.
+!   qui n'etait pas correcte.
 !
             IF (tnew.LT.RTT) THEN
@@ -4465 +5043 @@
       END
 
+!=====================================================================
+       SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas                                    &
+     &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
+     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
+     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                   &
+     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
+     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
+     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
+!
+     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
+     &         ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
+     &         ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                                          &
+     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
+     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
+     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
+ 
+       implicit none
+ 
+#include "dimensions.h"
+
+!-------------------------------------------------------------------------
+! Vertical interpolation of generic case forcing data onto mod_casel levels
+!-------------------------------------------------------------------------
+ 
+       integer nlevmax
+       parameter (nlevmax=41)
+       integer nlev_cas,mxcalc
+!       real play(llm), plev_prof(nlevmax) 
+!       real t_prof(nlevmax),q_prof(nlevmax)
+!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
+!       real ht_prof(nlevmax),vt_prof(nlevmax)
+!       real hq_prof(nlevmax),vq_prof(nlevmax)
+ 
+       real play(llm), plev_prof_cas(nlev_cas) 
+       real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
+       real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
+       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
+       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
+       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
+       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
+       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
+       real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
+       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
+ 
+       real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
+       real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
+       real u_mod_cas(llm),v_mod_cas(llm)
+       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
+       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
+       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
+       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
+       real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
+       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
+ 
+       integer l,k,k1,k2
+       real frac,frac1,frac2,fact
+ 
+       do l = 1, llm
+       print *,'debut interp2, play=',l,play(l)
+       enddo
+!      do l = 1, nlev_cas
+!      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
+!      enddo
+
+       do l = 1, llm
+
+        if (play(l).ge.plev_prof_cas(nlev_cas)) then
+ 
+        mxcalc=l
+        print *,'debut interp2, mxcalc=',mxcalc
+         k1=0
+         k2=0
+
+         if (play(l).le.plev_prof_cas(1)) then
+
+         do k = 1, nlev_cas-1
+          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
+            k1=k
+            k2=k+1
+          endif
+         enddo
+
+         if (k1.eq.0 .or. k2.eq.0) then
+          write(*,*) 'PB! k1, k2 = ',k1,k2
+          write(*,*) 'l,play(l) = ',l,play(l)/100
+         do k = 1, nlev_cas-1
+          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
+         enddo
+         endif
+
+         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
+         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
+         theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
+         thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
+         thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
+         qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
+         ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
+         qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
+         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
+         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
+         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
+         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
+         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
+         omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
+         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
+         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
+         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
+         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
+         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
+         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
+         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
+         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
+         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
+         dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
+         hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
+         vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
+         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
+         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
+         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
+     
+         else !play>plev_prof_cas(1)
+
+         k1=1
+         k2=2
+         print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
+         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
+         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
+         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
+         theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
+         thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
+         thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
+         qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
+         ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
+         qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
+         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
+         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
+         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
+         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
+         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
+         omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
+         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
+         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
+         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
+         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
+         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
+         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
+         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
+         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
+         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
+         dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
+         hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
+         vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
+         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
+         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
+         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
+
+         endif ! play.le.plev_prof_cas(1)
+
+        else ! above max altitude of forcing file
+ 
+!jyg
+         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
+         fact = max(fact,0.)                                           !jyg
+         fact = exp(-fact)                                             !jyg
+         t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
+         theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
+         thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
+         thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
+         qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
+         ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
+         qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
+         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
+         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
+         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                     !jyg
+         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                     !jyg
+         w_mod_cas(l)= 0.0                                             !jyg
+         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact 
+         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
+         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
+         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
+         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
+         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
+         dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
+         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
+         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
+         dth_mod_cas(l)= dth_prof_cas(nlev_cas) 
+         hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
+         vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
+         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact 
+         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
+         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
+ 
+        endif ! play
+ 
+       enddo ! l
+
+!       do l = 1,llm
+!       print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
+!     $        l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
+!       enddo
+ 
+          return
+          end
+!***************************************************************************** 
+
+