Index: /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3692)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3693)
@@ -190,6 +190,6 @@
       real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
       real :: d_u_nudge(llm),d_v_nudge(llm)
-      real :: du_adv(llm),dv_adv(llm)
-      real :: du_age(llm),dv_age(llm)
+!      real :: d_u_adv(llm),d_v_adv(llm)
+      real :: d_u_age(llm),d_v_age(llm)
       real :: alpha
       real :: ttt
@@ -273,8 +273,8 @@
       d_u_nudge(:)=0.
       d_v_nudge(:)=0.
-      du_adv(:)=0.
-      dv_adv(:)=0.
-      du_age(:)=0.
-      dv_age(:)=0.
+      d_u_adv(:)=0.
+      d_v_adv(:)=0.
+      d_u_age(:)=0.
+      d_v_age(:)=0.
       
 ! Initialization of Common turb_forcing
@@ -840,9 +840,4 @@
 
 #include "1D_interp_cases.h"
-   ! Vertical advection
-!  call lstendH(llm,nqtot,omega,d_t_vert_adv,d_q_vert_adv,q,temp,u,v,play)
-!  print*,'B d_t_adv ',d_t_adv(1:20)*86400
-!  print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400
-!  print*,'B dt omega ',omega
 
 !---------------------------------------------------------------------
@@ -872,9 +867,13 @@
       teta=temp*(pzero/play)**rkappa
       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_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 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))
       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)
@@ -938,90 +937,33 @@
        fcoriolis=2.*sin(rpi*xlat/180.)*romega
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!       if (forcing_radconv .or. forcing_fire) then
-!         fcoriolis=0.0
-!         dt_cooling=0.0
-!         d_t_adv=0.0
-!         d_q_adv=0.0
-!       endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!      if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
-!    &    .or.forcing_amma .or. forcing_type.eq.101) then
-!        fcoriolis=0.0 ; ug=0. ; vg=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!      if(forcing_rico) then
-!         dt_cooling=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!CRio:Attention modif sp??cifique cas de Caroline
-!     if (forcing_type==-1) then
-!        fcoriolis=0.
-!       
-!on calcule dt_cooling
-!       do l=1,llm
-!       if (play(l).ge.20000.) then
-!           dt_cooling(l)=-1.5/86400.
-!       elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
-!           dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
-!       else
-!           dt_cooling(l)=-1.*(temp(l)-200.)/86400.
-!       endif
-!       enddo
-!
-!     endif     
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!     if (forcing_sandu) then
-!        ug(1:llm)=u_mod(1:llm)
-!        vg(1:llm)=v_mod(1:llm)
-!     endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
       IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
                                    fcoriolis, xlat,mxcalc
 
-!       print *,'u-ug=',u-ug
-
-!!!!!!!!!!!!!!!!!!!!!!!!
-! Geostrophic wind
-! Le calcul ci dessous est insuffisamment precis
-!      du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
-!      dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
-!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Geostrophic forcing 
+!---------------------------------------------------------------------
+
+      IF ( forc_geo == 0 ) THEN
+              d_u_age(1:mxcalc)=0.
+              d_v_age(1:mxcalc)=0.
+      ELSE
        sfdt = sin(0.5*fcoriolis*timestep)
        cfdt = cos(0.5*fcoriolis*timestep)
 !       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
 !
-        du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
+        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
      &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
      &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
 !!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
 !
-       dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
+       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
      &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
      &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
 !!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
-!
-!!!!!!!!!!!!!!!!!!!!!!!!
+      ENDIF
+!
+!---------------------------------------------------------------------
 !  Nudging
-!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
       d_t_nudge(:) = 0.
       d_u_nudge(:) = 0.
@@ -1039,4 +981,7 @@
       ENDDO
 
+!---------------------------------------------------------------------
+!  Optional outputs
+!---------------------------------------------------------------------
 #ifdef OUTPUT_PHYS_SCM
       CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
@@ -1056,62 +1001,13 @@
 #endif
 
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!       if (forcing_fire) THEN
-!               print*,'Enlever cette section rapidement'
-!               stop
-!               
-!
-!!let ww=if ( alt le 1100 ) then alt*-0.00001 else 0
-!!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt)  else 0
-!!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt)  else 0
-!           d_t_adv=0.
-!           d_q_adv=0.
-!           teta=temp*(pzero/play)**rkappa
-!           d_t_adv=0.
-!           d_q_adv=0.
-!           do l=2,llm-1
-!              if (zlay(l)<=1100) then
-!                  wwww=-0.00001*zlay(l)
-!                  d_t_adv(l)=-wwww*(teta(l)-teta(l+1))/(zlay(l)-zlay(l+1)) /(pzero/play(l))**rkappa
-!                  d_q_adv(l,1:2)=-wwww*(q(l,1:2)-q(l+1,1:2))/(zlay(l)-zlay(l+1))
-!                  d_t_adv(l)=d_t_adv(l)+min(-3.75e-5 , -7.5e-8*zlay(l))
-!                  d_q_adv(l,1)=d_q_adv(l,1)+max( 1.5e-8 , 3e-11*zlay(l))
-!              endif
-!           enddo
-!
-!        endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!         call  writefield_phy('dv_age' ,dv_age,llm)
-!         call  writefield_phy('du_age' ,du_age,llm)
-!         call  writefield_phy('du_phys' ,du_phys,llm)
-!         call  writefield_phy('u_tend' ,u,llm)
-!         call  writefield_phy('u_g' ,ug,llm)
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! Increment state variables
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
-! au dessus de 700hpa, on relaxe vers les profils initiaux
-!     if (forcing_sandu .OR. forcing_astex) then
-!#include "1D_nudge_sandu_astex.h"
-!      else
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
      &              du_phys(1:mxcalc)                                       &
-     &             +du_age(1:mxcalc)+du_adv(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)                                       &
-     &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
+     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
      &             +d_v_nudge(1:mxcalc) )
         q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
@@ -1125,6 +1021,6 @@
      &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
            print* ,'dv_phys=',dv_phys
-           print* ,'dv_age=',dv_age
-           print* ,'dv_adv=',dv_adv
+           print* ,'d_v_age=',d_v_age
+           print* ,'d_v_adv=',d_v_adv
            print* ,'d_v_nudge=',d_v_nudge
            print*, v
