source: LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90 @ 2306

Last change on this file since 2306 was 2255, checked in by jyg, 10 years ago

Changes to pbl_surface and other routines concerning split/no-split.
+ pbl_surface_mod.F90: call cdrag for (w) region.
+ phyredem.F90: write wake_delta_pbl_TKE.
+ phys_output_write_mod.F90: control output of wake_delta_pbl_TKE by
IF(iflag_pbl_split>=1).
+ lmdz1d.F90: initialize wake_delta_pbl_TKE=0.
+ phys_output_ctrlout_mod.F90: suppression of accents in some variable
attributes.
+ cva_driver.F90: suppression of a print introduced in version 2253.

File size: 40.1 KB
Line 
1!#ifdef CPP_1D
2!#include "../dyn3d/mod_const_mpi.F90"
3!#include "../dyn3d_common/control_mod.F90"
4!#include "../dyn3d_common/infotrac.F90"
5!#include "../dyn3d_common/disvert.F90"
6
7
8      PROGRAM lmdz1d
9
10      USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
11      use phys_state_var_mod
12      use dimphy
13      use surface_data, only : type_ocean,ok_veget
14      use pbl_surface_mod, only : ftsoil, pbl_surface_init,                     &
15     &                            pbl_surface_final
16      use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
17
18      use infotrac ! new
19      use control_mod
20      USE indice_sol_mod
21      USE phyaqua_mod
22      USE mod_1D_cases_read
23      USE mod_1D_amma_read
24
25      implicit none
26#include "dimensions.h"
27#include "YOMCST.h"
28#include "temps.h"
29!!#include "control.h"
30#include "iniprint.h"
31#include "clesphys.h"
32#include "dimsoil.h"
33!#include "indicesol.h"
34
35#include "comvert.h"
36#include "compar1d.h"
37#include "flux_arp.h"
38#include "tsoilnudge.h"
39#include "fcg_gcssold.h"
40!!!#include "fbforcing.h"
41
42!=====================================================================
43! DECLARATIONS
44!=====================================================================
45
46!---------------------------------------------------------------------
47!  Externals
48!---------------------------------------------------------------------
49      external fq_sat
50      real fq_sat
51
52!---------------------------------------------------------------------
53!  Arguments d' initialisations de la physique (USER DEFINE)
54!---------------------------------------------------------------------
55
56      integer, parameter :: ngrid=1
57      real :: zcufi    = 1.
58      real :: zcvfi    = 1.
59
60!-      real :: nat_surf
61!-      logical :: ok_flux_surf
62!-      real :: fsens
63!-      real :: flat
64!-      real :: tsurf
65!-      real :: rugos
66!-      real :: qsol(1:2)
67!-      real :: qsurf
68!-      real :: psurf
69!-      real :: zsurf
70!-      real :: albedo
71!-
72!-      real :: time     = 0.
73!-      real :: time_ini
74!-      real :: xlat
75!-      real :: xlon
76!-      real :: wtsurf
77!-      real :: wqsurf
78!-      real :: restart_runoff
79!-      real :: xagesno
80!-      real :: qsolinp
81!-      real :: zpicinp
82!-
83      real :: fnday
84      real :: day, daytime
85      real :: day1
86      real :: heure
87      integer :: jour
88      integer :: mois
89      integer :: an
90 
91!---------------------------------------------------------------------
92!  Declarations related to forcing and initial profiles
93!---------------------------------------------------------------------
94
95        integer :: kmax = llm
96        integer llm700,nq1,nq2
97        INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
98        real timestep, frac
99        real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
100        real  uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
101        real  ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
102        real  dqtdxls(nlev_max),dqtdyls(nlev_max)
103        real  dqtdtls(nlev_max),thlpcar(nlev_max)
104        real  qprof(nlev_max,nqmx)
105
106!        integer :: forcing_type
107        logical :: forcing_les     = .false.
108        logical :: forcing_armcu   = .false.
109        logical :: forcing_rico    = .false.
110        logical :: forcing_radconv = .false.
111        logical :: forcing_toga    = .false.
112        logical :: forcing_twpice  = .false.
113        logical :: forcing_amma    = .false.
114        logical :: forcing_dice    = .false.
115        logical :: forcing_GCM2SCM = .false.
116        logical :: forcing_GCSSold = .false.
117        logical :: forcing_sandu   = .false.
118        logical :: forcing_astex   = .false.
119        logical :: forcing_fire    = .false.
120        logical :: forcing_case    = .false.
121        integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
122!                                                            (cf read_tsurf1d.F)
123
124!vertical advection computation
125!       real d_t_z(llm), d_q_z(llm)
126!       real d_t_dyn_z(llm), d_q_dyn_z(llm)
127!       real zz(llm)
128!       real zfact
129
130!flag forcings
131        logical :: nudge_wind=.true.
132        logical :: nudge_thermo=.false.
133        logical :: cptadvw=.true.
134!=====================================================================
135! DECLARATIONS FOR EACH CASE
136!=====================================================================
137!
138#include "1D_decl_cases.h"
139!
140!---------------------------------------------------------------------
141!  Declarations related to nudging
142!---------------------------------------------------------------------
143     integer :: nudge_max
144     parameter (nudge_max=9)
145     integer :: inudge_RHT=1
146     integer :: inudge_UV=2
147     logical :: nudge(nudge_max)
148     real :: t_targ(llm)
149     real :: rh_targ(llm)
150     real :: u_targ(llm)
151     real :: v_targ(llm)
152!
153!---------------------------------------------------------------------
154!  Declarations related to vertical discretization:
155!---------------------------------------------------------------------
156      real :: pzero=1.e5
157      real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
158      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
159
160!---------------------------------------------------------------------
161!  Declarations related to variables
162!---------------------------------------------------------------------
163
164      real :: phi(llm)
165      real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
166      real :: rlat_rad(1),rlon_rad(1)
167      real :: omega(llm+1),omega2(llm),rho(llm+1)
168      real :: ug(llm),vg(llm),fcoriolis
169      real :: sfdt, cfdt
170      real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
171      real :: dt_dyn(llm)
172      real :: dt_cooling(llm),d_th_adv(llm),d_t_nudge(llm)
173      real :: d_u_nudge(llm),d_v_nudge(llm)
174      real :: du_adv(llm),dv_adv(llm)
175      real :: alpha
176      real :: ttt
177
178      REAL, ALLOCATABLE, DIMENSION(:,:):: q
179      REAL, ALLOCATABLE, DIMENSION(:,:):: dq
180      REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
181      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
182      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
183!      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
184
185!---------------------------------------------------------------------
186!  Initialization of surface variables
187!---------------------------------------------------------------------
188      real :: run_off_lic_0(1)
189      real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
190      real :: tsoil(1,nsoilmx,nbsrf)
191!     real :: agesno(1,nbsrf)
192
193!---------------------------------------------------------------------
194!  Call to phyredem
195!---------------------------------------------------------------------
196      logical :: ok_writedem =.true.
197      real :: sollw_in = 0.
198      real :: solsw_in = 0.
199     
200!---------------------------------------------------------------------
201!  Call to physiq
202!---------------------------------------------------------------------
203      logical :: firstcall=.true.
204      logical :: lastcall=.false.
205      real :: phis    = 0.0
206      real :: dpsrf
207
208!---------------------------------------------------------------------
209!  Initializations of boundary conditions
210!---------------------------------------------------------------------
211      integer, parameter :: yd = 360
212      real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
213      real :: phy_alb (yd)  ! Albedo land only (old value condsurf_jyg=0.3)
214      real :: phy_sst (yd)  ! SST (will not be used; cf read_tsurf1d.F)
215      real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
216      real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
217      real :: phy_ice (yd) = 0.0 ! Fraction de glace
218      real :: phy_fter(yd) = 0.0 ! Fraction de terre
219      real :: phy_foce(yd) = 0.0 ! Fraction de ocean
220      real :: phy_fsic(yd) = 0.0 ! Fraction de glace
221      real :: phy_flic(yd) = 0.0 ! Fraction de glace
222
223!---------------------------------------------------------------------
224!  Fichiers et d'autres variables
225!---------------------------------------------------------------------
226      integer :: k,l,i,it=1,mxcalc
227      integer jcode
228      integer jjmp1
229      parameter (jjmp1=jjm+1-1/jjm)
230      REAL dudyn(iim+1,jjmp1,llm)
231      INTEGER read_climoz
232!Al1
233      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
234      data ecrit_slab_oc/-1/
235
236!=====================================================================
237! INITIALIZATIONS
238!=====================================================================
239! Initialization of Common turb_forcing
240       dtime_frcg = 0.
241       Turb_fcg_gcssold=.false.
242       hthturb_gcssold = 0.
243       hqturb_gcssold = 0.
244
245!---------------------------------------------------------------------
246! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
247!---------------------------------------------------------------------
248!Al1
249        call conf_unicol
250!Al1 moves this gcssold var from common fcg_gcssold to
251        Turb_fcg_gcssold = xTurb_fcg_gcssold
252! --------------------------------------------------------------------
253        close(1)
254!Al1
255        write(*,*) 'lmdz1d.def lu => unicol.def'
256
257! forcing_type defines the way the SCM is forced:
258!forcing_type = 0 ==> forcing_les = .true.
259!             initial profiles from file prof.inp.001
260!             no forcing by LS convergence ;
261!             surface temperature imposed ;
262!             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
263!forcing_type = 1 ==> forcing_radconv = .true.
264!             idem forcing_type = 0, but the imposed radiative cooling
265!             is set to 0 (hence, if iflag_radia=0 in physiq.def,
266!             then there is no radiative cooling at all)
267!forcing_type = 2 ==> forcing_toga = .true.
268!             initial profiles from TOGA-COARE IFA files
269!             LS convergence and SST imposed from TOGA-COARE IFA files
270!forcing_type = 3 ==> forcing_GCM2SCM = .true.
271!             initial profiles from the GCM output
272!             LS convergence imposed from the GCM output
273!forcing_type = 4 ==> forcing_twpice = .true.
274!             initial profiles from TWP-ICE cdf file
275!             LS convergence, omega and SST imposed from TWP-ICE files
276!forcing_type = 5 ==> forcing_rico = .true.
277!             initial profiles from RICO files
278!             LS convergence imposed from RICO files
279!forcing_type = 6 ==> forcing_amma = .true.
280!             initial profiles from AMMA nc file
281!             LS convergence, omega and surface fluxes imposed from AMMA file 
282!forcing_type = 7 ==> forcing_dice = .true.
283!             initial profiles and large scale forcings in dice_driver.nc
284!             Different stages: soil model alone, atm. model alone
285!             then both models coupled
286!forcing_type = 10 ==> forcing_case = .true.
287!             initial profiles and large scale forcings in cas.nc
288!             LS convergence, omega and SST imposed from CINDY-DYNAMO files
289!forcing_type = 40 ==> forcing_GCSSold = .true.
290!             initial profile from GCSS file
291!             LS convergence imposed from GCSS file
292!forcing_type = 50 ==> forcing_fire = .true.
293!             forcing from fire.nc
294!forcing_type = 59 ==> forcing_sandu = .true.
295!             initial profiles from sanduref file: see prof.inp.001
296!             SST varying with time and divergence constante: see ifa_sanduref.txt file
297!             Radiation has to be computed interactively
298!forcing_type = 60 ==> forcing_astex = .true.
299!             initial profiles from file: see prof.inp.001
300!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
301!             Radiation has to be computed interactively
302!forcing_type = 61 ==> forcing_armcu = .true.
303!             initial profiles from file: see prof.inp.001
304!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
305!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
306!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
307!             Radiation to be switched off
308!
309      if (forcing_type .eq.0) THEN
310       forcing_les = .true.
311      elseif (forcing_type .eq.1) THEN
312       forcing_radconv = .true.
313      elseif (forcing_type .eq.2) THEN
314       forcing_toga    = .true.
315      elseif (forcing_type .eq.3) THEN
316       forcing_GCM2SCM = .true.
317      elseif (forcing_type .eq.4) THEN
318       forcing_twpice = .true.
319      elseif (forcing_type .eq.5) THEN
320       forcing_rico = .true.
321      elseif (forcing_type .eq.6) THEN
322       forcing_amma = .true.
323      elseif (forcing_type .eq.7) THEN
324       forcing_dice = .true.
325      elseif (forcing_type .eq.10) THEN
326       forcing_case = .true.
327      elseif (forcing_type .eq.40) THEN
328       forcing_GCSSold = .true.
329      elseif (forcing_type .eq.50) THEN
330       forcing_fire = .true.
331      elseif (forcing_type .eq.59) THEN
332       forcing_sandu   = .true.
333      elseif (forcing_type .eq.60) THEN
334       forcing_astex   = .true.
335      elseif (forcing_type .eq.61) THEN
336       forcing_armcu = .true.
337       IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
338      else
339       write (*,*) 'ERROR : unknown forcing_type ', forcing_type
340       stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
341      ENDIF
342      print*,"forcing type=",forcing_type
343
344! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
345! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
346! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
347! through the common sst_forcing.
348
349        type_ts_forcing = 0
350        if (forcing_toga.or.forcing_sandu.or.forcing_astex .or. forcing_dice)                 &
351     &    type_ts_forcing = 1
352!
353! Initialization of the logical switch for nudging
354     jcode = iflag_nudge
355     do i = 1,nudge_max
356       nudge(i) = mod(jcode,10) .ge. 1
357       jcode = jcode/10
358     enddo
359!---------------------------------------------------------------------
360!  Definition of the run
361!---------------------------------------------------------------------
362
363      call conf_gcm( 99, .TRUE. )
364!-----------------------------------------------------------------------
365!   Choix du calendrier
366!   -------------------
367
368!      calend = 'earth_365d'
369      if (calend == 'earth_360d') then
370        call ioconf_calendar('360d')
371        write(*,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
372      else if (calend == 'earth_365d') then
373        call ioconf_calendar('noleap')
374        write(*,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
375      else if (calend == 'earth_366d') then
376        call ioconf_calendar('all_leap')
377        write(*,*)'CALENDRIER CHOISI: Terrestre bissextile'
378      else if (calend == 'gregorian') then
379        call ioconf_calendar('gregorian') ! not to be used by normal users
380        write(*,*)'CALENDRIER CHOISI: Gregorien'
381      else
382        write (*,*) 'ERROR : unknown calendar ', calend
383        stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
384      endif
385!-----------------------------------------------------------------------
386!
387!c Date :
388!      La date est supposee donnee sous la forme [annee, numero du jour dans
389!      l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
390!      On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
391!      Le numero du jour est dans "day". L heure est traitee separement.
392!      La date complete est dans "daytime" (l'unite est le jour).
393      if (nday>0) then
394         fnday=nday
395      else
396         fnday=-nday/float(day_step)
397      endif
398
399! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
400      IF(forcing_type .EQ. 61) fnday=53100./86400.
401! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
402      IF(forcing_type .EQ. 6) fnday=64800./86400.
403!     IF(forcing_type .EQ. 6) fnday=50400./86400.
404      annee_ref = anneeref
405      mois = 1
406      day_ref = dayref
407      heure = 0.
408      itau_dyn = 0
409      itau_phy = 0
410      call ymds2ju(annee_ref,mois,day_ref,heure,day)
411      day_ini = int(day)
412      day_end = day_ini + fnday
413
414      IF (forcing_type .eq.2) THEN
415! Convert the initial date of Toga-Coare to Julian day
416      call ymds2ju                                                          &
417     & (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
418
419      ELSEIF (forcing_type .eq.4) THEN
420! Convert the initial date of TWPICE to Julian day
421      call ymds2ju                                                          &
422     & (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi              &
423     & ,day_ju_ini_twpi)
424      ELSEIF (forcing_type .eq.6) THEN
425! Convert the initial date of AMMA to Julian day
426      call ymds2ju                                                          &
427     & (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma              &
428     & ,day_ju_ini_amma)
429      ELSEIF (forcing_type .eq.7) THEN
430! Convert the initial date of DICE to Julian day
431      call ymds2ju                                                         &
432     & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             &
433     & ,day_ju_ini_dice)
434      ELSEIF (forcing_type .eq.10) THEN
435! Convert the initial date to Julian day
436      print*,'time cindy',year_ini_cas,mth_ini_cas,day_ini_cas
437      call ymds2ju                                                         &
438     & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas &
439     & ,day_ju_ini_cas)
440      print*,'time cindy 2',day_ju_ini_cas
441      ELSEIF (forcing_type .eq.59) THEN
442! Convert the initial date of Sandu case to Julian day
443      call ymds2ju                                                          &
444     &   (year_ini_sandu,mth_ini_sandu,day_ini_sandu,                       &
445     &    time_ini*3600.,day_ju_ini_sandu)
446
447      ELSEIF (forcing_type .eq.60) THEN
448! Convert the initial date of Astex case to Julian day
449      call ymds2ju                                                          &
450     &   (year_ini_astex,mth_ini_astex,day_ini_astex,                        &
451     &    time_ini*3600.,day_ju_ini_astex)
452
453      ELSEIF (forcing_type .eq.61) THEN
454! Convert the initial date of Arm_cu case to Julian day
455      call ymds2ju                                                          &
456     & (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu          &
457     & ,day_ju_ini_armcu)
458      ENDIF
459
460      daytime = day + time_ini/24. ! 1st day and initial time of the simulation
461! Print out the actual date of the beginning of the simulation :
462      call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
463      print *,' Time of beginning : ',                                      &
464     &        year_print, month_print, day_print, sec_print
465
466!---------------------------------------------------------------------
467! Initialization of dimensions, geometry and initial state
468!---------------------------------------------------------------------
469      call init_phys_lmdz(1,1,llm,1,(/1/))
470      call suphel
471      call infotrac_init
472
473      if (nqtot>nqmx) STOP'Augmenter nqmx dans lmdz1d.F'
474      allocate(q(llm,nqtot)) ; q(:,:)=0.
475      allocate(dq(llm,nqtot))
476      allocate(dq_dyn(llm,nqtot))
477      allocate(d_q_adv(llm,nqtot))
478      allocate(d_q_nudge(llm,nqtot))
479!      allocate(d_th_adv(llm))
480
481!
482!   No ozone climatology need be read in this pre-initialization
483!          (phys_state_var_init is called again in physiq)
484      read_climoz = 0
485!
486      call phys_state_var_init(read_climoz)
487
488      if (ngrid.ne.klon) then
489         print*,'stop in inifis'
490         print*,'Probleme de dimensions :'
491         print*,'ngrid = ',ngrid
492         print*,'klon  = ',klon
493         stop
494      endif
495!!!=====================================================================
496!!! Feedback forcing values for Gateaux differentiation (al1)
497!!!=====================================================================
498!!! Surface Planck forcing bracketing call radiation
499!!      surf_Planck = 0.
500!!      surf_Conv   = 0.
501!!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
502!!! a mettre dans le lmdz1d.def ou autre
503!!
504!!
505      qsol = qsolinp
506      qsurf = fq_sat(tsurf,psurf/100.)
507      rlat=xlat
508      rlon=xlon
509      day1= day_ini
510      time=daytime-day
511      ts_toga(1)=tsurf ! needed by read_tsurf1d.F
512      rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf))
513
514!
515!! mpl et jyg le 22/08/2012 :
516!!  pour que les cas a flux de surface imposes marchent
517      IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
518       fsens=-wtsurf*rcpd*rho(1)
519       flat=-wqsurf*rlvtt*rho(1)
520       print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
521      ENDIF
522      print*,'Flux sol ',fsens,flat
523!!      ok_flux_surf=.false.
524!!      fsens=-wtsurf*rcpd*rho(1)
525!!      flat=-wqsurf*rlvtt*rho(1)
526!!!!
527
528! Vertical discretization and pressure levels at half and mid levels:
529
530      pa   = 5e4
531!!      preff= 1.01325e5
532      preff = psurf
533      IF (ok_old_disvert) THEN
534        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
535        print *,'On utilise disvert0'
536      ELSE
537        call disvert()
538        print *,'On utilise disvert'
539!       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
540!       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
541      ENDIF
542      sig_s=presnivs/preff
543      plev =ap+bp*psurf
544      play = 0.5*(plev(1:llm)+plev(2:llm+1))
545!cc      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
546
547      IF (forcing_type .eq. 59) THEN
548! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
549      write(*,*) '***********************'
550      do l = 1, llm
551       write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
552       if (trouve_700 .and. play(l).le.70000) then
553         llm700=l
554         print *,'llm700,play=',llm700,play(l)/100.
555         trouve_700= .false.
556       endif
557      enddo
558      write(*,*) '***********************'
559      ENDIF
560
561!
562!=====================================================================
563! EVENTUALLY, READ FORCING DATA :
564!=====================================================================
565
566#include "1D_read_forc_cases.h"
567
568      if (forcing_GCM2SCM) then
569        write (*,*) 'forcing_GCM2SCM not yet implemented'
570        stop 'in initialization'
571      endif ! forcing_GCM2SCM
572
573      print*,'mxcalc=',mxcalc
574      print*,'zlay=',zlay(mxcalc)
575      print*,'play=',play(mxcalc)
576
577!Al1 pour SST forced, appellé depuis ocean_forced_noice
578      ts_cur = tsurf ! SST used in read_tsurf1d
579!=====================================================================
580! Initialisation de la physique :
581!=====================================================================
582
583!  Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
584!
585! day_step, iphysiq lus dans gcm.def ci-dessus
586! timestep: calcule ci-dessous from rday et day_step
587! ngrid=1
588! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
589! rday: defini dans suphel.F (86400.)
590! day_ini: lu dans run.def (dayref)
591! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
592! airefi,zcufi,zcvfi initialises au debut de ce programme
593! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
594      day_step = float(nsplit_phys)*day_step/float(iphysiq)
595      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
596      timestep =rday/day_step
597      dtime_frcg = timestep
598!
599      zcufi=airefi
600      zcvfi=airefi
601!
602      rlat_rad(:)=rlat(:)*rpi/180.
603      rlon_rad(:)=rlon(:)*rpi/180.
604
605      call iniphysiq(iim,jjm,llm,rday,day_ini,timestep,                        &
606     &     rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
607      print*,'apres iniphysiq'
608
609! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
610      co2_ppm= 330.0
611      solaire=1370.0
612
613! Ecriture du startphy avant le premier appel a la physique.
614! On le met juste avant pour avoir acces a tous les champs
615
616      if (ok_writedem) then
617
618!--------------------------------------------------------------------------
619! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
620! need : qsol fder snow qsurf evap rugos agesno ftsoil
621!--------------------------------------------------------------------------
622
623        type_ocean = "force"
624        run_off_lic_0(1) = restart_runoff
625        call fonte_neige_init(run_off_lic_0)
626
627        fder=0.
628        snsrf(1,:)=0.        ! couverture de neige des sous surface
629        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
630        fevap=0.
631        z0m(1,:)=rugos     ! couverture de neige des sous surface
632        z0h(1,:)=rugos     ! couverture de neige des sous surface
633        agesno  = xagesno
634        tsoil(:,:,:)=tsurf
635!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
636!       tsoil(1,1,1)=299.18
637!       tsoil(1,2,1)=300.08
638!       tsoil(1,3,1)=301.88
639!       tsoil(1,4,1)=305.48
640!       tsoil(1,5,1)=308.00
641!       tsoil(1,6,1)=308.00
642!       tsoil(1,7,1)=308.00
643!       tsoil(1,8,1)=308.00
644!       tsoil(1,9,1)=308.00
645!       tsoil(1,10,1)=308.00
646!       tsoil(1,11,1)=308.00
647!-----------------------------------------------------------------------
648        call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
649
650!------------------ prepare limit conditions for limit.nc -----------------
651!--   Ocean force
652
653        print*,'avant phyredem'
654        pctsrf(1,:)=0.
655        if (nat_surf.eq.0.) then
656          pctsrf(1,is_oce)=1.
657          pctsrf(1,is_ter)=0.
658        else
659          pctsrf(1,is_oce)=0.
660          pctsrf(1,is_ter)=1.
661        end if
662
663        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
664     &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
665
666        zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
667        zpic = zpicinp
668        ftsol=tsurf
669        nsw=6 ! on met le nb de bandes SW=6, pour initialiser
670              ! 6 albedo, mais on peut quand meme tourner avec
671              ! moins. Seules les 2 ou 4 premiers seront lus
672        falb_dir=albedo
673        falb_dif=albedo
674        rugoro=rugos
675        t_ancien(1,:)=temp(:)
676        q_ancien(1,:)=q(:,1)
677        pbl_tke(:,:,:)=1.e-8
678        wake_delta_pbl_tke(:,:,:)=0.
679
680        rain_fall=0.
681        snow_fall=0.
682        solsw=0.
683        sollw=0.
684        radsol=0.
685        rnebcon=0.
686        ratqs=0.
687        clwcon=0.
688        zmea=0.
689        zstd=0.
690        zsig=0.
691        zgam=0.
692        zval=0.
693        zthe=0.
694        sig1=0.
695        w01=0.
696        u_ancien(1,:)=u(:)
697        v_ancien(1,:)=v(:)
698 
699!------------------------------------------------------------------------
700! Make file containing restart for the physics (startphy.nc)
701!
702! NB: List of the variables to be written by phyredem (via put_field):
703! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
704! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
705! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
706! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
707! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
708! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
709! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
710! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
711!------------------------------------------------------------------------
712!Al1 =============== restart option ==========================
713        if (.not.restart) then
714          call phyredem ("startphy.nc")
715        else
716! (desallocations)
717        print*,'callin surf final'
718          call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
719        print*,'after surf final'
720          CALL fonte_neige_final(run_off_lic_0)
721        endif
722
723        ok_writedem=.false.
724        print*,'apres phyredem'
725
726      endif ! ok_writedem
727     
728!------------------------------------------------------------------------
729! Make file containing boundary conditions (limit.nc) **Al1->restartdyn***
730! --------------------------------------------------
731! NB: List of the variables to be written in limit.nc
732!     (by writelim.F, subroutine of 1DUTILS.h):
733!        phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
734!        phy_fter,phy_foce,phy_flic,phy_fsic)
735!------------------------------------------------------------------------
736      do i=1,yd
737        phy_nat(i)  = nat_surf
738        phy_alb(i)  = albedo
739        phy_sst(i)  = tsurf ! read_tsurf1d will be used instead
740        phy_rug(i)  = rugos
741        phy_fter(i) = pctsrf(1,is_ter)
742        phy_foce(i) = pctsrf(1,is_oce)
743        phy_fsic(i) = pctsrf(1,is_sic)
744        phy_flic(i) = pctsrf(1,is_lic)
745      enddo
746
747! fabrication de limit.nc
748      call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,             &
749     &               phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
750
751
752      call phys_state_var_end
753!Al1
754      if (restart) then
755        print*,'call to restart dyn 1d'
756        Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs,          &
757     &              u,v,temp,q,omega2)
758
759       print*,'fnday,annee_ref,day_ref,day_ini',                            &
760     &     fnday,annee_ref,day_ref,day_ini
761!**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
762       day = day_ini
763       day_end = day_ini + nday
764       daytime = day + time_ini/24. ! 1st day and initial time of the simulation
765
766! Print out the actual date of the beginning of the simulation :
767       call ju2ymds(daytime, an, mois, jour, heure)
768       print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
769
770       day = int(daytime)
771       time=daytime-day
772 
773       print*,'****** intialised fields from restart1dyn *******'
774       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
775       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
776       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
777! raz for safety
778       do l=1,llm
779         dq_dyn(l,1) = 0.
780       enddo
781      endif
782!Al1 ================  end restart =================================
783      IF (ecrit_slab_oc.eq.1) then
784         open(97,file='div_slab.dat',STATUS='UNKNOWN')
785       elseif (ecrit_slab_oc.eq.0) then
786         open(97,file='div_slab.dat',STATUS='OLD')
787       endif
788!
789!---------------------------------------------------------------------
790!    Initialize target profile for RHT nudging if needed
791!---------------------------------------------------------------------
792      if (nudge(inudge_RHT)) then
793        call nudge_RHT_init(plev,play,temp,q(:,1),t_targ,rh_targ)
794      endif
795      if (nudge(inudge_UV)) then
796        call nudge_UV_init(plev,play,u,v,u_targ,v_targ)
797      endif
798!
799!=====================================================================
800! START OF THE TEMPORAL LOOP :
801!=====================================================================
802           
803      do while(it.le.nint(fnday*day_step))
804
805       if (prt_level.ge.1) then
806         print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',                       &
807     &                                it,day,time,nint(fnday*day_step)
808         print*,'PAS DE TEMPS ',timestep
809       endif
810!Al1 demande de restartphy.nc
811       if (it.eq.nint(fnday*day_step)) lastcall=.True.
812
813!---------------------------------------------------------------------
814! Interpolation of forcings in time and onto model levels
815!---------------------------------------------------------------------
816
817#include "1D_interp_cases.h"
818
819      if (forcing_GCM2SCM) then
820        write (*,*) 'forcing_GCM2SCM not yet implemented'
821        stop 'in time loop'
822      endif ! forcing_GCM2SCM
823
824!---------------------------------------------------------------------
825!  Geopotential :
826!---------------------------------------------------------------------
827
828        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
829        do l = 1, llm-1
830          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
831     &    (play(l)-play(l+1))/(play(l)+play(l+1))
832        enddo
833
834!---------------------------------------------------------------------
835! Listing output for debug prt_level>=1
836!---------------------------------------------------------------------
837       if (prt_level>=1) then
838         print *,' avant physiq : -------- day time ',day,time
839         write(*,*) 'firstcall,lastcall,phis',                               &
840     &               firstcall,lastcall,phis
841         write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',                   &
842     &        'presniv','plev','play','phi'
843         write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,                   &
844     &         presnivs(l),plev(l),play(l),phi(l),l=1,llm)
845         write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',                    &
846     &         'presniv','u','v','temp','q1','q2','omega2'
847         write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,         &
848     &   presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
849       endif
850
851!---------------------------------------------------------------------
852!   Call physiq :
853!---------------------------------------------------------------------
854
855        call physiq(ngrid,llm,                                              &
856     &              firstcall,lastcall,                                     &
857     &              day,time,timestep,                                      &
858     &              plev,play,phi,phis,presnivs,                            &
859     &              u,v,temp,q,omega2,                                      &
860     &              du_phys,dv_phys,dt_phys,dq,dpsrf,                        &
861     &              dudyn)
862        firstcall=.false.
863
864!---------------------------------------------------------------------
865! Listing output for debug prt_level>=1
866!---------------------------------------------------------------------
867        if (prt_level>=1) then
868          write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',                  &
869     &        'presniv','plev','play','phi'
870          write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,                  &
871     &    presnivs(l),plev(l),play(l),phi(l),l=1,llm)
872          write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',                   &
873     &         'presniv','u','v','temp','q1','q2','omega2'
874          write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,       &
875     &    presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
876          write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',                   &
877     &         'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'   
878           write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,            &
879     &      presnivs(l),86400*du_phys(l),86400*dv_phys(l),                   &
880     &       86400*dt_phys(l),86400*dq(l,1),dq(l,2),l=1,llm)
881          write(*,*) 'dpsrf',dpsrf
882        endif
883!---------------------------------------------------------------------
884!   Add physical tendencies :
885!---------------------------------------------------------------------
886
887       fcoriolis=2.*sin(rpi*xlat/180.)*romega
888       if (forcing_radconv .or. forcing_fire) then
889         fcoriolis=0.0
890         dt_cooling=0.0
891         d_th_adv=0.0
892         d_q_adv=0.0
893       endif
894
895       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
896     &    .or.forcing_amma) then
897         fcoriolis=0.0 ; ug=0. ; vg=0.
898       endif
899         if(forcing_rico) then
900          dt_cooling=0.
901        endif
902
903      print*, 'fcoriolis ', fcoriolis, xlat
904
905        du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
906       dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
907
908!!!!!!!!!!!!!!!!!!!!!!!!
909! Geostrophic wind
910!!!!!!!!!!!!!!!!!!!!!!!!
911       sfdt = sin(0.5*fcoriolis*timestep)
912       cfdt = cos(0.5*fcoriolis*timestep)
913!
914        du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
915     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
916     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
917!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
918!
919       dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
920     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
921     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
922!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
923!
924!!!!!!!!!!!!!!!!!!!!!!!!
925!  Nudging
926!!!!!!!!!!!!!!!!!!!!!!!!
927      d_t_nudge(:) = 0.
928      d_q_nudge(:,:) = 0.
929      d_u_nudge(:) = 0.
930      d_v_nudge(:) = 0.
931      if (nudge(inudge_RHT)) then
932        call nudge_RHT(timestep,plev,play,t_targ,rh_targ,temp,q(:,1),     &
933    &                  d_t_nudge,d_q_nudge(:,1))
934      endif
935      if (nudge(inudge_UV)) then
936        call nudge_UV(timestep,plev,play,u_targ,v_targ,u,v,     &
937    &                  d_u_nudge,d_v_nudge)
938      endif
939!
940!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
941!         call  writefield_phy('dv_age' ,dv_age,llm)
942!         call  writefield_phy('du_age' ,du_age,llm)
943!         call  writefield_phy('du_phys' ,du_phys,llm)
944!         call  writefield_phy('u_tend' ,u,llm)
945!         call  writefield_phy('u_g' ,ug,llm)
946!
947!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
948!! Increment state variables
949!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
950! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
951! au dessus de 700hpa, on relaxe vers les profils initiaux
952      if (forcing_sandu .OR. forcing_astex) then
953#include "1D_nudge_sandu_astex.h"
954      else
955        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
956     &              du_phys(1:mxcalc)                                       &
957     &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
958     &             +d_u_nudge(1:mxcalc) )           
959        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
960     &              dv_phys(1:mxcalc)                                       &
961     &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
962     &             +d_v_nudge(1:mxcalc) )
963        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
964     &                dq(1:mxcalc,:)                                        &
965     &               +d_q_adv(1:mxcalc,:)                                   &
966     &               +d_q_nudge(1:mxcalc,:) )
967
968        if (prt_level.ge.1) then
969          print *,                                                          &
970     &    'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',         &
971     &              temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
972           print*,du_phys
973           print*, v
974           print*, vg
975        endif
976
977        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
978     &              dt_phys(1:mxcalc)                                       &
979     &             +d_th_adv(1:mxcalc)                                      &
980     &             +d_t_nudge(1:mxcalc)                                      &
981     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
982
983      endif  ! forcing_sandu or forcing_astex
984
985        teta=temp*(pzero/play)**rkappa
986!
987!---------------------------------------------------------------------
988!   Nudge soil temperature if requested
989!---------------------------------------------------------------------
990
991      IF (nudge_tsoil) THEN
992       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
993     &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
994      ENDIF
995
996!---------------------------------------------------------------------
997!   Add large-scale tendencies (advection, etc) :
998!---------------------------------------------------------------------
999
1000!cc nrlmd
1001!cc        tmpvar=teta
1002!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1003!cc
1004!cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
1005!cc        tmpvar(:)=q(:,1)
1006!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1007!cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
1008!cc        tmpvar(:)=q(:,2)
1009!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1010!cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
1011
1012!---------------------------------------------------------------------
1013!   Air temperature :
1014!---------------------------------------------------------------------       
1015        if (lastcall) then
1016          print*,'Pas de temps final ',it
1017          call ju2ymds(daytime, an, mois, jour, heure)
1018          print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1019        endif
1020
1021!  incremente day time
1022!        print*,'daytime bef',daytime,1./day_step
1023        daytime = daytime+1./day_step
1024!Al1dbg
1025        day = int(daytime+0.1/day_step)
1026!        time = max(daytime-day,0.0)
1027!Al1&jyg: correction de bug
1028!cc        time = real(mod(it,day_step))/day_step
1029        time = time_ini/24.+real(mod(it,day_step))/day_step
1030!        print*,'daytime nxt time',daytime,time
1031        it=it+1
1032
1033      enddo
1034
1035!Al1
1036      if (ecrit_slab_oc.ne.-1) close(97)
1037
1038!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1039! -------------------------------------
1040       call dyn1dredem("restart1dyn.nc",                                    &
1041     &              plev,play,phi,phis,presnivs,                            &
1042     &              u,v,temp,q,omega2)
1043
1044        CALL abort_gcm ('lmdz1d   ','The End  ',0)
1045
1046      end
1047
1048#include "1DUTILS.h"
1049#include "1Dconv.h"
1050
1051!#endif
1052
Note: See TracBrowser for help on using the repository browser.