Ignore:
Timestamp:
Oct 22, 2020, 2:50:18 PM (4 years ago)
Author:
evignon
Message:

Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90

    r3693 r3780  
    7575      real :: zcufi    = 1.
    7676      real :: zcvfi    = 1.
    77 
    78 !-      real :: nat_surf
    79 !-      logical :: ok_flux_surf
    80 !-      real :: fsens
    81 !-      real :: flat
    82 !-      real :: tsurf
    83 !-      real :: rugos
    84 !-      real :: qsol(1:2)
    85 !-      real :: qsurf
    86 !-      real :: psurf
    87 !-      real :: zsurf
    88 !-      real :: albedo
    89 !-
    90 !-      real :: time     = 0.
    91 !-      real :: time_ini
    92 !-      real :: xlat
    93 !-      real :: xlon
    94 !-      real :: wtsurf
    95 !-      real :: wqsurf
    96 !-      real :: restart_runoff
    97 !-      real :: xagesno
    98 !-      real :: qsolinp
    99 !-      real :: zpicinp
    100 !-
    10177      real :: fnday
    10278      real :: day, daytime
    … …  
    141117        logical :: forcing_case2   = .false.
    142118        logical :: forcing_SCM   = .false.
    143         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
    144 !                                                            (cf read_tsurf1d.F)
    145119
    146120!flag forcings
    … …  
    148122        logical :: nudge_thermo=.false.
    149123        logical :: cptadvw=.true.
     124
     125
    150126!=====================================================================
    151127! DECLARATIONS FOR EACH CASE
    … …  
    248224!
    249225      integer :: it_end ! iteration number of the last call
    250 !Al1
     226!Al1,plev,play,phi,phis,presnivs,
    251227      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
    252228      data ecrit_slab_oc/-1/
    … …  
    278254      d_v_age(:)=0.
    279255     
     256     
    280257! Initialization of Common turb_forcing
    281258       dtime_frcg = 0.
    … …  
    290267! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
    291268!---------------------------------------------------------------------
    292 !Al1
    293269        call conf_unicol
    294270!Al1 moves this gcssold var from common fcg_gcssold to
    … …  
    296272! --------------------------------------------------------------------
    297273        close(1)
    298 !Al1
    299274        write(*,*) 'lmdz1d.def lu => unicol.def'
    300275
    … …  
    302277       year_ini_cas=1997
    303278       ! It is possible that those parameters are run twice.
    304 
    305279       ! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
     280
     281
    306282       call getin('anneeref',year_ini_cas)
    307283       call getin('dayref',day_deb)
    … …  
    309285       call getin('time_ini',heure_ini_cas)
    310286
    311         type_ts_forcing = 0
    312         IF (nat_surf==0) type_ts_forcing=1 ! SST forcee sur OCEAN
    313         print*,'NATURE DE LA SURFACE ',nat_surf
     287      print*,'NATURE DE LA SURFACE ',nat_surf
    314288!
    315289! Initialization of the logical switch for nudging
     290
    316291     jcode = iflag_nudge
    317292     do i = 1,nudge_max
    … …  
    319294       jcode = jcode/10
    320295     enddo
    321 !---------------------------------------------------------------------
     296!-----------------------------------------------------------------------
    322297!  Definition of the run
    323 !---------------------------------------------------------------------
     298!-----------------------------------------------------------------------
    324299
    325300      call conf_gcm( 99, .TRUE. )
    … …  
    343318      allocate( phy_flic(year_len)) ! Fraction de glace
    344319      phy_flic(:)=0.0
     320
     321
    345322!-----------------------------------------------------------------------
    346323!   Choix du calendrier
    … …  
    373350!      Le numero du jour est dans "day". L heure est traitee separement.
    374351!      La date complete est dans "daytime" (l'unite est le jour).
     352
     353
    375354      if (nday>0) then
    376355         fnday=nday
    … …  
    409388! Initialization of dimensions, geometry and initial state
    410389!---------------------------------------------------------------------
    411 !      call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
     390!     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
    412391!     but we still need to initialize dimphy module (klon,klev,etc.)  here.
    413392      call init_dimphy1D(1,llm)
    … …  
    446425!!! Feedback forcing values for Gateaux differentiation (al1)
    447426!!!=====================================================================
    448 !!! Surface Planck forcing bracketing call radiation
    449 !!      surf_Planck = 0.
    450 !!      surf_Conv   = 0.
    451 !!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
    452 !!! a mettre dans le lmdz1d.def ou autre
    453 !!
    454427!!
    455428      qsol = qsolinp
    … …  
    469442      ENDIF
    470443      print*,'Flux sol ',fsens,flat
    471 !!      ok_flux_surf=.false.
    472 !!      fsens=-wtsurf*rcpd*rho(1)
    473 !!      flat=-wqsurf*rlvtt*rho(1)
    474 !!!!
    475444
    476445! Vertical discretization and pressure levels at half and mid levels:
    … …  
    496465      plev =ap+bp*psurf
    497466      play = 0.5*(plev(1:llm)+plev(2:llm+1))
    498       zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
     467      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles.
    499468
    500469      IF (forcing_type .eq. 59) THEN
    … …  
    527496      print*,'mxcalc=',mxcalc
    528497!     print*,'zlay=',zlay(mxcalc)
    529       print*,'play=',play(mxcalc)
    530 
    531 !Al1 pour SST forced, appell?? depuis ocean_forced_noice
    532       ts_cur = tsurf ! SST used in read_tsurf1d
     498!      print*,'play=',play(mxcalc)
     499
     500!! When surface temperature is forced
     501      tg= tsurf ! surface T used in read_tsurf1d
     502
     503
    533504!=====================================================================
    534505! Initialisation de la physique :
    … …  
    546517! airefi,zcufi,zcvfi initialises au debut de ce programme
    547518! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
     519
     520
    548521      day_step = float(nsplit_phys)*day_step/float(iphysiq)
    549522      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
    … …  
    563536     ! e.g. for cell boundaries, which are meaningless in 1D; so pad these
    564537     ! with '0.' when necessary
     538
    565539      call iniphysiq(iim,jjm,llm, &
    566540           1,comm_lmdz, &
    … …  
    653627              ! 6 albedo, mais on peut quand meme tourner avec
    654628              ! moins. Seules les 2 ou 4 premiers seront lus
     629              print*, 'les albedos sont', albedo
    655630        falb_dir=albedo
    656631        falb_dif=albedo
    … …  
    664639        prw_ancien = 0.
    665640!jyg<
    666 !!        pbl_tke(:,:,:)=1.e-8
     641!!      pbl_tke(:,:,:)=1.e-8
    667642        pbl_tke(:,:,:)=0.
    668         pbl_tke(:,2,:)=1.e-2
     643! EV: pourquoi????
     644!        pbl_tke(:,2,:)=1.e-2
    669645        PRINT *, ' pbl_tke dans lmdz1d '
    670646        if (prt_level .ge. 5) then
    … …  
    675651
    676652!>jyg
    677 
    678653        rain_fall=0.
    679654        snow_fall=0.
    … …  
    715690        v_ancien(1,:)=v(:)
    716691 
    717 u10m=0.
    718 v10m=0.
    719 ale_wake=0.
    720 ale_bl_stat=0.
     692        u10m=0.
     693        v10m=0.
     694        ale_wake=0.
     695        ale_bl_stat=0.
    721696
    722697!------------------------------------------------------------------------
    … …  
    738713! to be set at some arbitratry convenient values.
    739714!------------------------------------------------------------------------
    740 !Al1 =============== restart option ==========================
     715!Al1 =============== restart option ======================================
    741716        if (.not.restart) then
    742717          iflag_pbl = 5
    … …  
    803778       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
    804779       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
    805        print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
     780       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
    806781! raz for safety
    807782       do l=1,llm
    … …  
    809784       enddo
    810785      endif
    811 !Al1 ================  end restart =================================
     786!======================  end restart =================================
    812787      IF (ecrit_slab_oc.eq.1) then
    813788         open(97,file='div_slab.dat',STATUS='UNKNOWN')
    … …  
    820795       CALL iophys_ini
    821796#endif
     797
     798!=====================================================================
    822799! START OF THE TEMPORAL LOOP :
    823800!=====================================================================
    824801           
    825802      it_end = nint(fnday*day_step)
    826 !test JLD     it_end = 10
    827803      do while(it.le.it_end)
    828804
    … …  
    832808         print*,'PAS DE TEMPS ',timestep
    833809       endif
    834 !Al1 demande de restartphy.nc
    835810       if (it.eq.it_end) lastcall=.True.
    836811
    … …  
    844819!  Geopotential :
    845820!---------------------------------------------------------------------
    846 
     821!        phis(1)=zsurf*RG
     822!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
    847823        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
     824
    848825        do l = 1, llm-1
    849826          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
    850827     &    (play(l)-play(l+1))/(play(l)+play(l+1))
    851828        enddo
     829
    852830
    853831!---------------------------------------------------------------------
    … …  
    950928       sfdt = sin(0.5*fcoriolis*timestep)
    951929       cfdt = cos(0.5*fcoriolis*timestep)
    952 !       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
    953 !
     930
    954931        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
    955932     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
    … …  
    10301007        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
    10311008     &              dt_phys(1:mxcalc)                                       &
    1032      &             +d_t_adv(1:mxcalc)                                      &
    1033      &             +d_t_nudge(1:mxcalc)                                      &
     1009     &             +d_t_adv(1:mxcalc)                                       &
     1010     &             +d_t_nudge(1:mxcalc)                                     &
    10341011     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
    10351012
    10361013
    1037 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1014!=======================================================================
    10381015!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    1039 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1040 !     endif  ! forcing_sandu or forcing_astex
    1041 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1016!=======================================================================
    10421017
    10431018        teta=temp*(pzero/play)**rkappa
    1044 !
     1019
    10451020!---------------------------------------------------------------------
    10461021!   Nudge soil temperature if requested
    … …  
    10801055
    10811056!  incremente day time
    1082 !        print*,'daytime bef',daytime,1./day_step
    10831057        daytime = daytime+1./day_step
    1084 !Al1dbg
    10851058        day = int(daytime+0.1/day_step)
    10861059!        time = max(daytime-day,0.0)
    … …  
    10881061!cc        time = real(mod(it,day_step))/day_step
    10891062        time = time_ini/24.+real(mod(it,day_step))/day_step
    1090 !        print*,'daytime nxt time',daytime,time
    10911063        it=it+1
    10921064
    10931065      enddo
    10941066
    1095 !Al1
    10961067      if (ecrit_slab_oc.ne.-1) close(97)
    10971068
    10981069!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
    1099 ! -------------------------------------
     1070! ---------------------------------------------------------------------------
    11001071       call dyn1dredem("restart1dyn.nc",                                    &
    11011072     &              plev,play,phi,phis,presnivs,                            &
Note: See TracChangeset for help on using the changeset viewer.