Changeset 3851 for LMDZ6/branches/LMDZ-tracers/libf/phylmd/dyn1d/scm.F90

Timestamp:

Feb 22, 2021, 12:44:07 PM (4 years ago)

Author:

dcugnet

Message:

Update the branch to the current trunk.

Location:

LMDZ6/branches/LMDZ-tracers

Files:

• . (modified) (1 prop)
• libf/phylmd/dyn1d/scm.F90 (modified) (33 diffs)

LMDZ6/branches/LMDZ-tracers/libf/phylmd/dyn1d/scm.F90

@@ -75 +75 @@
       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 +117 @@
         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 +122 @@
         logical :: nudge_thermo=.false.
         logical :: cptadvw=.true.
+
+
 !=====================================================================
 ! DECLARATIONS FOR EACH CASE
@@ -248 +224 @@
 !
       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 +254 @@
       d_v_age(:)=0.
       
+      
 ! Initialization of Common turb_forcing
        dtime_frcg = 0.
@@ -290 +267 @@
 ! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
 !---------------------------------------------------------------------
-!Al1
         call conf_unicol
 !Al1 moves this gcssold var from common fcg_gcssold to 
@@ -296 +272 @@
 ! --------------------------------------------------------------------
         close(1)
-!Al1
         write(*,*) 'lmdz1d.def lu => unicol.def'
 
@@ -302 +277 @@
        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 +285 @@
        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 +294 @@
        jcode = jcode/10
      enddo
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 !  Definition of the run
-!---------------------------------------------------------------------
+!-----------------------------------------------------------------------
 
       call conf_gcm( 99, .TRUE. )
@@ -343 +318 @@
       allocate( phy_flic(year_len)) ! Fraction de glace
       phy_flic(:)=0.0
+
+
 !-----------------------------------------------------------------------
 !   Choix du calendrier
@@ -373 +350 @@
 !      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 +388 @@
 ! 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)
@@ -433 +412 @@
 !          (phys_state_var_init is called again in physiq)
       read_climoz = 0
-!
+      nsw=6
+
       call phys_state_var_init(read_climoz)
 
@@ -446 +426 @@
 !!! 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 +443 @@
       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 +466 @@
       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 +497 @@
       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 +518 @@
 ! 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 +537 @@
      ! 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, &
@@ -650 +625 @@
         zpic = zpicinp
         ftsol=tsurf
-        nsw=6 ! on met le nb de bandes SW=6, pour initialiser
-              ! 6 albedo, mais on peut quand meme tourner avec
-              ! moins. Seules les 2 ou 4 premiers seront lus
         falb_dir=albedo
         falb_dif=albedo
@@ -664 +636 @@
         prw_ancien = 0.
 !jyg<
-!!        pbl_tke(:,:,:)=1.e-8
-        pbl_tke(:,:,:)=0.
-        pbl_tke(:,2,:)=1.e-2
-        PRINT *, ' pbl_tke dans lmdz1d '
-        if (prt_level .ge. 5) then 
-         DO nsrf = 1,4
-           PRINT *,'pbl_tke(1,:,',nsrf,') ',pbl_tke(1,:,nsrf)
-         ENDDO
-        end if
-
+! 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.
@@ -715 +680 @@
         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 +703 @@
 ! to be set at some arbitratry convenient values.
 !------------------------------------------------------------------------
-!Al1 =============== restart option ==========================
+!Al1 =============== restart option ======================================
         if (.not.restart) then
           iflag_pbl = 5
@@ -803 +768 @@
        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 +774 @@
        enddo
       endif
-!Al1 ================  end restart =================================
+!======================  end restart =================================
       IF (ecrit_slab_oc.eq.1) then
          open(97,file='div_slab.dat',STATUS='UNKNOWN')
@@ -820 +785 @@
        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 +798 @@
          print*,'PAS DE TEMPS ',timestep
        endif
-!Al1 demande de restartphy.nc
        if (it.eq.it_end) lastcall=.True.
 
@@ -844 +809 @@
 !  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 +918 @@
        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 +997 @@
         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 +1045 @@
 
 !  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 +1051 @@
 !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,                            &