Index: LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3693)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3780)
@@ -75,28 +75,4 @@
       real :: zcufi    = 1.
       real :: zcvfi    = 1.
-
-!-      real :: nat_surf
-!-      logical :: ok_flux_surf
-!-      real :: fsens
-!-      real :: flat
-!-      real :: tsurf
-!-      real :: rugos
-!-      real :: qsol(1:2)
-!-      real :: qsurf
-!-      real :: psurf
-!-      real :: zsurf
-!-      real :: albedo
-!-
-!-      real :: time     = 0.
-!-      real :: time_ini
-!-      real :: xlat 
-!-      real :: xlon 
-!-      real :: wtsurf 
-!-      real :: wqsurf 
-!-      real :: restart_runoff 
-!-      real :: xagesno 
-!-      real :: qsolinp 
-!-      real :: zpicinp 
-!-
       real :: fnday 
       real :: day, daytime 
@@ -141,6 +117,4 @@
         logical :: forcing_case2   = .false.
         logical :: forcing_SCM   = .false.
-        integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
-!                                                            (cf read_tsurf1d.F)
 
 !flag forcings
@@ -148,4 +122,6 @@
         logical :: nudge_thermo=.false.
         logical :: cptadvw=.true.
+
+
 !=====================================================================
 ! DECLARATIONS FOR EACH CASE
@@ -248,5 +224,5 @@
 !
       integer :: it_end ! iteration number of the last call
-!Al1
+!Al1,plev,play,phi,phis,presnivs,
       integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
       data ecrit_slab_oc/-1/
@@ -278,4 +254,5 @@
       d_v_age(:)=0.
       
+      
 ! Initialization of Common turb_forcing
        dtime_frcg = 0.
@@ -290,5 +267,4 @@
 ! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
 !---------------------------------------------------------------------
-!Al1
         call conf_unicol
 !Al1 moves this gcssold var from common fcg_gcssold to 
@@ -296,5 +272,4 @@
 ! --------------------------------------------------------------------
         close(1)
-!Al1
         write(*,*) 'lmdz1d.def lu => unicol.def'
 
@@ -302,6 +277,7 @@
        year_ini_cas=1997
        ! It is possible that those parameters are run twice.
-
        ! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
+
+
        call getin('anneeref',year_ini_cas)
        call getin('dayref',day_deb)
@@ -309,9 +285,8 @@
        call getin('time_ini',heure_ini_cas)
 
-        type_ts_forcing = 0
-        IF (nat_surf==0) type_ts_forcing=1 ! SST forcee sur OCEAN
-        print*,'NATURE DE LA SURFACE ',nat_surf
+      print*,'NATURE DE LA SURFACE ',nat_surf
 !
 ! Initialization of the logical switch for nudging
+
      jcode = iflag_nudge
      do i = 1,nudge_max
@@ -319,7 +294,7 @@
        jcode = jcode/10
      enddo
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 !  Definition of the run
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 
       call conf_gcm( 99, .TRUE. )
@@ -343,4 +318,6 @@
       allocate( phy_flic(year_len)) ! Fraction de glace
       phy_flic(:)=0.0
+
+
 !-----------------------------------------------------------------------
 !   Choix du calendrier
@@ -373,4 +350,6 @@
 !      Le numero du jour est dans "day". L heure est traitee separement.
 !      La date complete est dans "daytime" (l'unite est le jour).
+
+
       if (nday>0) then
          fnday=nday
@@ -409,5 +388,5 @@
 ! Initialization of dimensions, geometry and initial state
 !---------------------------------------------------------------------
-!      call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
+!     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
 !     but we still need to initialize dimphy module (klon,klev,etc.)  here.
       call init_dimphy1D(1,llm)
@@ -446,10 +425,4 @@
 !!! Feedback forcing values for Gateaux differentiation (al1)
 !!!=====================================================================
-!!! Surface Planck forcing bracketing call radiation
-!!      surf_Planck = 0.
-!!      surf_Conv   = 0.
-!!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
-!!! a mettre dans le lmdz1d.def ou autre
-!!
 !!
       qsol = qsolinp
@@ -469,8 +442,4 @@
       ENDIF
       print*,'Flux sol ',fsens,flat
-!!      ok_flux_surf=.false.
-!!      fsens=-wtsurf*rcpd*rho(1)
-!!      flat=-wqsurf*rlvtt*rho(1)
-!!!!
 
 ! Vertical discretization and pressure levels at half and mid levels:
@@ -496,5 +465,5 @@
       plev =ap+bp*psurf
       play = 0.5*(plev(1:llm)+plev(2:llm+1))
-      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
+      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles. 
 
       IF (forcing_type .eq. 59) THEN
@@ -527,8 +496,10 @@
       print*,'mxcalc=',mxcalc
 !     print*,'zlay=',zlay(mxcalc)
-      print*,'play=',play(mxcalc)
-
-!Al1 pour SST forced, appell?? depuis ocean_forced_noice
-      ts_cur = tsurf ! SST used in read_tsurf1d
+!      print*,'play=',play(mxcalc)
+
+!! When surface temperature is forced
+      tg= tsurf ! surface T used in read_tsurf1d
+
+
 !=====================================================================
 ! Initialisation de la physique : 
@@ -546,4 +517,6 @@
 ! airefi,zcufi,zcvfi initialises au debut de ce programme
 ! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
+
+
       day_step = float(nsplit_phys)*day_step/float(iphysiq)
       write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
@@ -563,4 +536,5 @@
      ! e.g. for cell boundaries, which are meaningless in 1D; so pad these 
      ! with '0.' when necessary
+
       call iniphysiq(iim,jjm,llm, &
            1,comm_lmdz, &
@@ -653,4 +627,5 @@
               ! 6 albedo, mais on peut quand meme tourner avec
               ! moins. Seules les 2 ou 4 premiers seront lus
+              print*, 'les albedos sont', albedo
         falb_dir=albedo
         falb_dif=albedo
@@ -664,7 +639,8 @@
         prw_ancien = 0.
 !jyg<
-!!        pbl_tke(:,:,:)=1.e-8
+!!      pbl_tke(:,:,:)=1.e-8
         pbl_tke(:,:,:)=0.
-        pbl_tke(:,2,:)=1.e-2
+! EV: pourquoi????
+!        pbl_tke(:,2,:)=1.e-2
         PRINT *, ' pbl_tke dans lmdz1d '
         if (prt_level .ge. 5) then 
@@ -675,5 +651,4 @@
 
 !>jyg
-
         rain_fall=0.
         snow_fall=0.
@@ -715,8 +690,8 @@
         v_ancien(1,:)=v(:)
  
-u10m=0.
-v10m=0.
-ale_wake=0.
-ale_bl_stat=0.
+        u10m=0.
+        v10m=0.
+        ale_wake=0.
+        ale_bl_stat=0.
 
 !------------------------------------------------------------------------
@@ -738,5 +713,5 @@
 ! to be set at some arbitratry convenient values.
 !------------------------------------------------------------------------
-!Al1 =============== restart option ==========================
+!Al1 =============== restart option ======================================
         if (.not.restart) then
           iflag_pbl = 5
@@ -803,5 +778,5 @@
        print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
        print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
-       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
+       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
 ! raz for safety
        do l=1,llm
@@ -809,5 +784,5 @@
        enddo
       endif
-!Al1 ================  end restart =================================
+!======================  end restart =================================
       IF (ecrit_slab_oc.eq.1) then
          open(97,file='div_slab.dat',STATUS='UNKNOWN')
@@ -820,9 +795,10 @@
        CALL iophys_ini
 #endif
+
+!=====================================================================
 ! START OF THE TEMPORAL LOOP :
 !=====================================================================
            
       it_end = nint(fnday*day_step)
-!test JLD     it_end = 10
       do while(it.le.it_end)
 
@@ -832,5 +808,4 @@
          print*,'PAS DE TEMPS ',timestep
        endif
-!Al1 demande de restartphy.nc
        if (it.eq.it_end) lastcall=.True.
 
@@ -844,10 +819,13 @@
 !  Geopotential :
 !---------------------------------------------------------------------
-
+!        phis(1)=zsurf*RG
+!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
         phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
+
         do l = 1, llm-1
           phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
      &    (play(l)-play(l+1))/(play(l)+play(l+1))
         enddo
+
 
 !---------------------------------------------------------------------
@@ -950,6 +928,5 @@
        sfdt = sin(0.5*fcoriolis*timestep)
        cfdt = cos(0.5*fcoriolis*timestep)
-!       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
-!
+
         d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
      &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
@@ -1030,17 +1007,15 @@
         temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
      &              dt_phys(1:mxcalc)                                       &
-     &             +d_t_adv(1:mxcalc)                                      &
-     &             +d_t_nudge(1:mxcalc)                                      &
+     &             +d_t_adv(1:mxcalc)                                       &
+     &             +d_t_nudge(1:mxcalc)                                     &
      &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!=======================================================================
 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!     endif  ! forcing_sandu or forcing_astex
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!=======================================================================
 
         teta=temp*(pzero/play)**rkappa
-!
+
 !---------------------------------------------------------------------
 !   Nudge soil temperature if requested
@@ -1080,7 +1055,5 @@
 
 !  incremente day time
-!        print*,'daytime bef',daytime,1./day_step
         daytime = daytime+1./day_step
-!Al1dbg
         day = int(daytime+0.1/day_step)
 !        time = max(daytime-day,0.0)
@@ -1088,14 +1061,12 @@
 !cc        time = real(mod(it,day_step))/day_step
         time = time_ini/24.+real(mod(it,day_step))/day_step
-!        print*,'daytime nxt time',daytime,time
         it=it+1
 
       enddo
 
-!Al1
       if (ecrit_slab_oc.ne.-1) close(97)
 
 !Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
-! -------------------------------------
+! ---------------------------------------------------------------------------
        call dyn1dredem("restart1dyn.nc",                                    &
      &              plev,play,phi,phis,presnivs,                            &
