Changeset 1795 for LMDZ5/branches/testing/libf/phy1d/lmdz1d.F

Timestamp:

Jul 18, 2013, 10:20:28 AM (11 years ago)

Author:

Ehouarn Millour

Message:

Version testing basee sur la r1794

---

Testing release based on r1794

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phy1d/lmdz1d.F (modified) (21 diffs)

LMDZ5/branches/testing/libf/phy1d/lmdz1d.F

@@ -6 +6 @@
       use dimphy
       use surface_data, only : type_ocean,ok_veget
-      use pbl_surface_mod, only : pbl_surface_init, pbl_surface_final
+      use pbl_surface_mod, only : ftsoil, pbl_surface_init,
+     $                            pbl_surface_final
       use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 
 
       use infotrac ! new
       use control_mod
+      USE indice_sol_mod
 
       implicit none
@@ -20 +22 @@
 #include "clesphys.h"
 #include "dimsoil.h"
-#include "indicesol.h"
+!#include "indicesol.h"
 
 #include "comvert.h"
 #include "compar1d.h"
 #include "flux_arp.h"
+#include "tsoilnudge.h"
 #include "fcg_gcssold.h"
 !!!#include "fbforcing.h"
@@ -86 +89 @@
 
         integer :: kmax = llm
-        integer nlev_max
-        parameter (nlev_max = 100)
+        integer nlev_max,llm700
+        parameter (nlev_max = 1000)
         real timestep, frac, timeit
         real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max),
@@ -98 +101 @@
 c        integer :: forcing_type
         logical :: forcing_les     = .false.
-        logical :: forcing_armcu  = .false.
+        logical :: forcing_armcu   = .false.
         logical :: forcing_rico    = .false.
         logical :: forcing_radconv = .false.
         logical :: forcing_toga    = .false.
         logical :: forcing_twpice  = .false.
+        logical :: forcing_amma    = .false.
         logical :: forcing_GCM2SCM = .false.
         logical :: forcing_GCSSold = .false.
+        logical :: forcing_sandu   = .false.
+        logical :: forcing_astex   = .false.
+        logical :: forcing_fire    = .false.
         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
 !                                                            (cf read_tsurf1d.F)
 
 !vertical advection computation
-        real d_t_z(llm), d_q_z(llm)
-        real d_t_dyn_z(llm), d_q_dyn_z(llm)
-        real zz(llm)
-        real zfact
+!       real d_t_z(llm), d_q_z(llm)
+!       real d_t_dyn_z(llm), d_q_dyn_z(llm)
+!       real zz(llm)
+!       real zfact
 
 !flag forcings
@@ -129 +136 @@
       real :: pzero=1.e5
       real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
-      real :: playd(llm),zlayd(llm)
+      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
 
 !---------------------------------------------------------------------
@@ -137 +144 @@
       integer :: iq
       real :: phi(llm)
-      real :: teta(llm),temp(llm),u(llm),v(llm)
+      real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
+      real :: rlat_rad(1),rlon_rad(1)
       real :: omega(llm+1),omega2(llm),rho(llm+1)
       real :: ug(llm),vg(llm),fcoriolis
+      real :: sfdt, cfdt
       real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
       real :: du_dyn(llm),dv_dyn(llm),dt_dyn(llm)
@@ -194 +203 @@
 !  Fichiers et d'autres variables
 !---------------------------------------------------------------------
-      real ttt
+      real ttt,bow,q1
       integer :: ierr,k,l,i,it=1,mxcalc
       integer jjmp1
@@ -250 +259 @@
 !             initial profiles from RICO files 
 !             LS convergence imposed from RICO files 
+!forcing_type = 6 ==> forcing_amma = .true.
+!             initial profiles from AMMA nc file 
+!             LS convergence, omega and surface fluxes imposed from AMMA file  
 !forcing_type = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
 !             LS convergence imposed from GCSS file
+!forcing_type = 59 ==> forcing_sandu = .true.
+!             initial profiles from sanduref file: see prof.inp.001
+!             SST varying with time and divergence constante: see ifa_sanduref.txt file
+!             Radiation has to be computed interactively
+!forcing_type = 60 ==> forcing_astex = .true.
+!             initial profiles from file: see prof.inp.001 
+!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
+!             Radiation has to be computed interactively
 !forcing_type = 61 ==> forcing_armcu = .true.
-!             initial profile from arm_cu file
-!             LS convergence imposed from arm_cu file
+!             initial profiles from file: see prof.inp.001 
+!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
+!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
+!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 
+!             Radiation to be switched off
 !
       if (forcing_type .eq.0) THEN
@@ -269 +292 @@
       elseif (forcing_type .eq.5) THEN
        forcing_rico = .true.
+      elseif (forcing_type .eq.6) THEN
+       forcing_amma = .true.
       elseif (forcing_type .eq.40) THEN
        forcing_GCSSold = .true.
+      elseif (forcing_type .eq.59) THEN
+       forcing_sandu   = .true.
+      elseif (forcing_type .eq.60) THEN
+       forcing_astex   = .true.
       elseif (forcing_type .eq.61) THEN
        forcing_armcu = .true.
@@ -276 +305 @@
       else
        write (*,*) 'ERROR : unknown forcing_type ', forcing_type
-       stop 'Forcing_type should be 0,1,2,3 or 40'
+       stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
       ENDIF
       print*,"forcing type=",forcing_type
@@ -286 +315 @@
 
         type_ts_forcing = 0
-        if (forcing_toga) type_ts_forcing = 1
+        if (forcing_toga .or. forcing_sandu .or. forcing_astex) 
+     :    type_ts_forcing = 1
 
 !---------------------------------------------------------------------
@@ -325 +355 @@
 c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
       IF(forcing_type .EQ. 61) fnday=53100./86400.
+c Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
+      IF(forcing_type .EQ. 6) fnday=64800./86400.
       annee_ref = anneeref
       mois = 1
@@ -334 +366 @@
       day_ini = day
       day_end = day_ini + nday
+
+      IF (forcing_type .eq.2) THEN
 ! Convert the initial date of Toga-Coare to Julian day
       call ymds2ju
      $ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
 
+      ELSEIF (forcing_type .eq.4) THEN
 ! Convert the initial date of TWPICE to Julian day
       call ymds2ju
      $ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
      $ ,day_ju_ini_twpi)
-
-! Convert the initial date of Arm_cu to Julian day
+      ELSEIF (forcing_type .eq.6) THEN
+! Convert the initial date of AMMA to Julian day
+      call ymds2ju
+     $ (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma
+     $ ,day_ju_ini_amma)
+
+      ELSEIF (forcing_type .eq.59) THEN
+! Convert the initial date of Sandu case to Julian day
+      call ymds2ju
+     $   (year_ini_sandu,mth_ini_sandu,day_ini_sandu,
+     $    time_ini*3600.,day_ju_ini_sandu)
+
+      ELSEIF (forcing_type .eq.60) THEN
+! Convert the initial date of Astex case to Julian day
+      call ymds2ju
+     $   (year_ini_astex,mth_ini_astex,day_ini_astex,
+     $    time_ini*3600.,day_ju_ini_astex)
+
+      ELSEIF (forcing_type .eq.61) THEN
+
+! Convert the initial date of Arm_cu case to Julian day
       call ymds2ju
      $ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
      $ ,day_ju_ini_armcu)
+      ENDIF
 
       daytime = day + time_ini/24. ! 1st day and initial time of the simulation 
@@ -418 +473 @@
 !!      preff= 1.01325e5
       preff = psurf
-      call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+      IF (ok_old_disvert) THEN
+        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+        print *,'On utilise disvert0'
+      ELSE
+        call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     :                 scaleheight)
+        print *,'On utilise disvert'
+c       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
+c       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
+      ENDIF
       sig_s=presnivs/preff
       plev =ap+bp*psurf
@@ -424 +488 @@
 ccc      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
 
+      IF (forcing_type .eq. 59) THEN
+! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
       write(*,*) '***********************'
       do l = 1, llm
        write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
+       if (trouve_700 .and. play(l).le.70000) then
+         llm700=l
+         print *,'llm700,play=',llm700,play(l)/100.
+         trouve_700= .false.
+       endif
       enddo
       write(*,*) '***********************'
+      ENDIF
 
 c
@@ -460 +532 @@
 ! rday: defini dans suphel.F (86400.)
 ! day_ini: lu dans run.def (dayref)
-! rlat,rlon lus dans lmdz1d.def
+! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
 ! airefi,zcufi,zcvfi initialises au debut de ce programme
 ! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
@@ -470 +542 @@
       zcufi=airefi
       zcvfi=airefi
+!
+      rlat_rad(:)=rlat(:)*rpi/180.
+      rlon_rad(:)=rlon(:)*rpi/180.
 
       call iniphysiq(ngrid,llm,rday,day_ini,timestep,
-     .     rlat,rlon,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,1)
+     .     rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
       print*,'apres iniphysiq'
 
@@ -501 +576 @@
         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(qsol, fder, snsrf, qsurfsrf,
      &                                    evap, frugs, agesno, tsoil)
@@ -734 +822 @@
        endif
 
-       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice) then
+       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice
+     :    .or.forcing_amma) then
          fcoriolis=0.0 ; ug=0. ; vg=0.
        endif
@@ -748 +837 @@
      : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
 
+!!!!!!!!!!!!!!!!!!!!!!!!
+! Geostrophic wind
+!!!!!!!!!!!!!!!!!!!!!!!!
+       sfdt = sin(0.5*fcoriolis*timestep)
+       cfdt = cos(0.5*fcoriolis*timestep)
+!
+        du_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*
+     :          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +
+     :           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
+!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!         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
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         u(1:mxcalc)=u(1:mxcalc) + timestep*(
      :              du_phys(1:mxcalc)
@@ -773 +888 @@
 
         teta=temp*(pzero/play)**rkappa
+!
+!---------------------------------------------------------------------
+!   Nudge soil temperature if requested
+!---------------------------------------------------------------------
+
+      IF (nudge_tsoil) THEN
+       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) 
+     .  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
+      ENDIF
 
 !---------------------------------------------------------------------