source: LMDZ5/branches/LF-private/libf/phy1d/lmdz1d.F @ 3839

Last change on this file since 3839 was 1866, checked in by Laurent Fairhead, 11 years ago

La distinction avec ou sans writelim n'a plus lieu d'être
MP Lefebvre


the writelim/nowritelim distinction is of no further use
MP Lefebvre

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