Index: LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90	(revision 5302)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90	(revision 5433)
@@ -1,8 +1,8 @@
 SUBROUTINE scm
 
-USE flux_arp_mod_h
-      USE compbl_mod_h
-         USE clesphys_mod_h
-      USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
+   USE flux_arp_mod_h
+   USE compbl_mod_h
+   USE clesphys_mod_h
+   USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
        clwcon, detr_therm, &
@@ -33,7 +33,5 @@
    USE indice_sol_mod
    USE phyaqua_mod
-!  USE mod_1D_cases_read
    USE mod_1D_cases_read_std
-   !USE mod_1D_amma_read
    USE print_control_mod, ONLY: lunout, prt_level
    USE iniphysiq_mod, ONLY: iniphysiq
@@ -49,9 +47,9 @@
    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
    USE dimsoil_mod_h, ONLY: nsoilmx
-  USE yomcst_mod_h
-  USE tsoilnudge_mod_h
-  USE fcg_gcssold_mod_h
-  USE compar1d_mod_h
-  USE date_cas_mod_h
+   USE yomcst_mod_h
+   USE tsoilnudge_mod_h
+   USE fcg_gcssold_mod_h
+   USE compar1d_mod_h
+   USE date_cas_mod_h
 implicit none
 
@@ -86,5 +84,5 @@
 
         integer :: kmax = llm
-        integer llm700,nq1,nq2
+        integer llm700,nq1,nq2,iflag_1d_vert_adv
         INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
         real timestep, frac
@@ -164,5 +162,4 @@
       real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
       real :: d_u_nudge(llm),d_v_nudge(llm)
-!      real :: d_u_adv(llm),d_v_adv(llm)
       real :: d_u_age(llm),d_v_age(llm)
       real :: alpha
@@ -174,5 +171,4 @@
       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
-!      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
 
 !---------------------------------------------------------------------
@@ -182,5 +178,4 @@
       real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
       real :: tsoil(1,nsoilmx,nbsrf)
-!     real :: agesno(1,nbsrf)
 
 !---------------------------------------------------------------------
@@ -220,10 +215,8 @@
       integer jcode
       INTEGER read_climoz 
-!
       integer :: it_end ! iteration number of the last call
-!Al1,plev,play,phi,phis,presnivs,
       integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
       data ecrit_slab_oc/-1/
-!
+
 !     if flag_inhib_forcing = 0, tendencies of forcing are added
 !                           <> 0, tendencies of forcing are not added
@@ -282,7 +275,6 @@
        mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
        call getin('time_ini',heure_ini_cas)
-
       print*,'NATURE DE LA SURFACE ',nat_surf
-!
+
 ! Initialization of the logical switch for nudging
 
@@ -292,4 +284,9 @@
        jcode = jcode/10
      enddo
+! numerical scheme for vertical advection
+
+     iflag_1d_vert_adv=1
+     call getin('iflag_1d_vert_adv',iflag_1d_vert_adv)
+
 !-----------------------------------------------------------------------
 !  Definition of the run
@@ -398,5 +395,4 @@
       allocate(d_q_adv(llm,nqtot)) 
       allocate(d_q_nudge(llm,nqtot)) 
-!      allocate(d_th_adv(llm)) 
 
       q(:,:) = 0.
@@ -406,5 +402,4 @@
       d_q_nudge(:,:) = 0.
 
-!
 !   No ozone climatology need be read in this pre-initialization
 !          (phys_state_var_init is called again in physiq)
@@ -424,5 +419,5 @@
 !!! Feedback forcing values for Gateaux differentiation (al1)
 !!!=====================================================================
-!!
+      
       qsol = qsolinp
       qsurf = fq_sat(tsurf,psurf/100.)
@@ -433,5 +428,4 @@
       rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf)) 
 
-!
 !! mpl et jyg le 22/08/2012 : 
 !!  pour que les cas a flux de surface imposes marchent
@@ -487,5 +481,4 @@
 
       INCLUDE "1D_read_forc_cases.h"
-   print*,'A d_t_adv ',d_t_adv(1:20)*86400
 
       if (forcing_GCM2SCM) then
@@ -567,16 +560,4 @@
         agesno  = xagesno
         tsoil(:,:,:)=tsurf
-!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
-!       tsoil(1,1,1)=299.18
-!       tsoil(1,2,1)=300.08
-!       tsoil(1,3,1)=301.88
-!       tsoil(1,4,1)=305.48
-!       tsoil(1,5,1)=308.00
-!       tsoil(1,6,1)=308.00
-!       tsoil(1,7,1)=308.00
-!       tsoil(1,8,1)=308.00
-!       tsoil(1,9,1)=308.00
-!       tsoil(1,10,1)=308.00
-!       tsoil(1,11,1)=308.00
 !-----------------------------------------------------------------------
         call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
@@ -635,10 +616,4 @@
           rvc_ancien = 0.
         ENDIF
-!jyg<
-! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
-!!      pbl_tke(:,:,:)=1.e-8
-!        pbl_tke(:,:,:)=0.
-!        pbl_tke(:,2,:)=1.e-2
-!>jyg
         rain_fall=0.
         snow_fall=0.
@@ -754,5 +729,4 @@
 
       call phys_state_var_end
-!Al1
       if (restart) then
         print*,'call to restart dyn 1d'
@@ -762,5 +736,4 @@
        print*,'fnday,annee_ref,day_ref,day_ini',                            &
      &     fnday,annee_ref,day_ref,day_ini
-!**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
        day = day_ini
        day_end = day_ini + nday
@@ -780,5 +753,5 @@
 ! raz for safety
        do l=1,llm
-         d_q_vert_adv(l,1) = 0.
+         d_q_vert_adv(l,:) = 0.
        enddo
       endif
@@ -819,5 +792,4 @@
 !---------------------------------------------------------------------
         phis(1)=zsurf*RG
-!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
 
         ! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
@@ -842,18 +814,48 @@
          z_adv=play
       ENDIF
-
-      teta=temp*(pzero/play)**rkappa
+      
+      ! vertical tendencies computed as d X / d t = -W  d X / d z
+
+      IF (iflag_1d_vert_adv .EQ. 0) THEN
+
+      ! old centered numerical scheme
       do l=2,llm-1
-        ! vertical tendencies computed as d X / d t = -W  d X / d z
         d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
         d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
         ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
-        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
-        d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))
+        d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l-1))/(z_adv(l+1)-z_adv(l-1)) 
+        d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1))
       enddo
-      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
-      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
-      d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
-      d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
+
+
+      ELSE IF (iflag_1d_vert_adv .EQ. 1) THEN
+      ! upstream numerical scheme
+
+      do l=2,llm-1
+        IF ( ( ( forc_w .EQ. 1 ) .AND. ( w_adv(l) .GT. 0. ) ) .OR. &
+             ( ( forc_w .NE. 1 ) .AND. ( w_adv(l) .LT. 0. ) ) ) THEN
+          d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1))
+          d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1))
+          d_t_vert_adv(l)=-w_adv(l)*(temp(l)-temp(l-1))/(z_adv(l)-z_adv(l-1))
+          d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1))
+        ELSE
+          d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l))
+          d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l))
+          d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l))/(z_adv(l+1)-z_adv(l))
+          d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l))
+        ENDIF
+      enddo
+
+
+      ENDIF ! numerical scheme for vertical advection
+      
+      
+      IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
+          u(1:mxcalc)=u(1:mxcalc) + timestep * d_u_vert_adv(1:mxcalc)
+          v(1:mxcalc)=v(1:mxcalc) + timestep * d_v_vert_adv(1:mxcalc)
+          q(1:mxcalc,:)=q(1:mxcalc,:) + timestep * d_q_vert_adv(1:mxcalc,:)
+          temp(1:mxcalc)=temp(1:mxcalc) + timestep * d_t_vert_adv(1:mxcalc)
+          teta=temp*(pzero/play)**rkappa
+      ENDIF
    
    ENDIF
@@ -940,6 +942,4 @@
 !---------------------------------------------------------------------
 !  Nudging
-!  EV: rewrite the section to account for a time- and height-varying
-!  nudging
 !---------------------------------------------------------------------
       d_t_nudge(:) = 0.
@@ -1000,4 +1000,30 @@
 
       ENDDO
+
+!-----------------------------------------------------------
+! horizontal forcings (advection) and nudging
+!-----------------------------------------------------------
+
+    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
+
+        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
+     &              du_phys(1:mxcalc)                                       &
+     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
+     &             +d_u_nudge(1:mxcalc) )            
+        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
+     &              dv_phys(1:mxcalc)                                       &
+     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
+     &             +d_v_nudge(1:mxcalc) )
+        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
+     &                dq(1:mxcalc,:)                                        &
+     &               +d_q_adv(1:mxcalc,:)                                   &
+     &               +d_q_nudge(1:mxcalc,:) )
+
+
+        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
+     &              dt_phys(1:mxcalc)                                       &
+     &             +d_t_adv(1:mxcalc)                                       &
+     &             +d_t_nudge(1:mxcalc)                                     &
+     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
 
 !---------------------------------------------------------------------
@@ -1022,41 +1048,7 @@
 END IF
 
-    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
-
-        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
-     &              du_phys(1:mxcalc)                                       &
-     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
-     &             +d_u_nudge(1:mxcalc) )            
-        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
-     &              dv_phys(1:mxcalc)                                       &
-     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
-     &             +d_v_nudge(1:mxcalc) )
-        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
-     &                dq(1:mxcalc,:)                                        &
-     &               +d_q_adv(1:mxcalc,:)                                   &
-     &               +d_q_nudge(1:mxcalc,:) )
-
-        if (prt_level.ge.3) then
-          print *,                                                          &
-     &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
-     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
-           print* ,'dv_phys=',dv_phys
-           print* ,'d_v_age=',d_v_age
-           print* ,'d_v_adv=',d_v_adv
-           print* ,'d_v_nudge=',d_v_nudge
-           print*, v
-           print*, vg
-        endif
-
-        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
-     &              dt_phys(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 !!
-!=======================================================================
+
+
+! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
 
         teta=temp*(pzero/play)**rkappa
@@ -1071,21 +1063,6 @@
       ENDIF
 
-!---------------------------------------------------------------------
-!   Add large-scale tendencies (advection, etc) :
-!---------------------------------------------------------------------
-
-!cc nrlmd
-!cc        tmpvar=teta
-!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
-!cc
-!cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
-!cc        tmpvar(:)=q(:,1)
-!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
-!cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
-!cc        tmpvar(:)=q(:,2)
-!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
-!cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
-
-   END IF ! end if tendency of tendency should be added
+
+   END IF ! end if in case tendency should be added
 
 !---------------------------------------------------------------------
@@ -1101,7 +1078,4 @@
         daytime = daytime+1./day_step
         day = int(daytime+0.1/day_step)
-!        time = max(daytime-day,0.0)
-!Al1&jyg: correction de bug
-!cc        time = real(mod(it,day_step))/day_step
         time = time_ini/24.+real(mod(it,day_step))/day_step
         it=it+1
