source: LMDZ5/trunk/libf/phy1d/lmdz1d.F @ 1969

Last change on this file since 1969 was 1969, checked in by fhourdin, 11 years ago

Correction. Declarations manquantes.
Bug fixing (mixing declarations).

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