source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/lmdz1d.F90

Last change on this file was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

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