source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/dyn1d/scm.F90 @ 5411

Last change on this file since 5411 was 4727, checked in by idelkadi, 14 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

File size: 40.5 KB
Line 
1SUBROUTINE scm
2
3   USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
4   USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
5       clwcon, detr_therm, &
6       qsol, fevap, z0m, z0h, agesno, &
7       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
8       falb_dir, falb_dif, &
9       ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
10       rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
11       solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, rneb_ancien, &
12       wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
13       wake_deltaq, wake_deltat, wake_s, wake_dens, &
14       awake_dens, cv_gen, wake_cstar, &
15       zgam, zmax0, zmea, zpic, zsig, &
16       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, ql_ancien, qs_ancien, &
17       prlw_ancien, prsw_ancien, prw_ancien, &
18       u10m,v10m,ale_wake,ale_bl_stat
19
20 
21   USE dimphy
22   USE surface_data, only : type_ocean,ok_veget
23   USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
24                                 pbl_surface_final
25   USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
26
27   USE infotrac ! new
28   USE control_mod
29   USE indice_sol_mod
30   USE phyaqua_mod
31!  USE mod_1D_cases_read
32   USE mod_1D_cases_read_std
33   !USE mod_1D_amma_read
34   USE print_control_mod, ONLY: lunout, prt_level
35   USE iniphysiq_mod, ONLY: iniphysiq
36   USE mod_const_mpi, ONLY: comm_lmdz
37   USE physiq_mod, ONLY: physiq
38   USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
39                          preff, aps, bps, pseudoalt, scaleheight
40   USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
41                        itau_dyn, itau_phy, start_time, year_len
42   USE phys_cal_mod, ONLY : year_len_phys_cal_mod => year_len
43
44      implicit none
45      INCLUDE "dimensions.h"
46      INCLUDE "YOMCST.h"
47!!      INCLUDE "control.h"
48      INCLUDE "clesphys.h"
49      INCLUDE "dimsoil.h"
50!      INCLUDE "indicesol.h"
51
52      INCLUDE "compar1d.h"
53      INCLUDE "flux_arp.h"
54      INCLUDE "date_cas.h"
55      INCLUDE "tsoilnudge.h"
56      INCLUDE "fcg_gcssold.h"
57      INCLUDE "compbl.h"
58
59!=====================================================================
60! DECLARATIONS
61!=====================================================================
62
63#undef OUTPUT_PHYS_SCM
64
65!---------------------------------------------------------------------
66!  Externals
67!---------------------------------------------------------------------
68      external fq_sat
69      real fq_sat
70
71!---------------------------------------------------------------------
72!  Arguments d' initialisations de la physique (USER DEFINE)
73!---------------------------------------------------------------------
74
75      integer, parameter :: ngrid=1
76      real :: zcufi    = 1.
77      real :: zcvfi    = 1.
78      real :: fnday
79      real :: day, daytime
80      real :: day1
81      real :: heure
82      integer :: jour
83      integer :: mois
84      integer :: an
85 
86!---------------------------------------------------------------------
87!  Declarations related to forcing and initial profiles
88!---------------------------------------------------------------------
89
90        integer :: kmax = llm
91        integer llm700,nq1,nq2
92        INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
93        real timestep, frac
94        real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
95        real  uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
96        real  ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
97        real  dqtdxls(nlev_max),dqtdyls(nlev_max)
98        real  dqtdtls(nlev_max),thlpcar(nlev_max)
99        real  qprof(nlev_max,nqmx)
100
101!        integer :: forcing_type
102        logical :: forcing_les     = .false.
103        logical :: forcing_armcu   = .false.
104        logical :: forcing_rico    = .false.
105        logical :: forcing_radconv = .false.
106        logical :: forcing_toga    = .false.
107        logical :: forcing_twpice  = .false.
108        logical :: forcing_amma    = .false.
109        logical :: forcing_dice    = .false.
110        logical :: forcing_gabls4  = .false.
111
112        logical :: forcing_GCM2SCM = .false.
113        logical :: forcing_GCSSold = .false.
114        logical :: forcing_sandu   = .false.
115        logical :: forcing_astex   = .false.
116        logical :: forcing_fire    = .false.
117        logical :: forcing_case    = .false.
118        logical :: forcing_case2   = .false.
119        logical :: forcing_SCM   = .false.
120
121!flag forcings
122        logical :: nudge_wind=.true.
123        logical :: nudge_thermo=.false.
124        logical :: cptadvw=.true.
125
126
127!=====================================================================
128! DECLARATIONS FOR EACH CASE
129!=====================================================================
130!
131      INCLUDE "1D_decl_cases.h"
132!
133!---------------------------------------------------------------------
134!  Declarations related to nudging
135!---------------------------------------------------------------------
136     integer :: nudge_max
137     parameter (nudge_max=9)
138     integer :: inudge_RHT=1
139     integer :: inudge_UV=2
140     logical :: nudge(nudge_max)
141     real :: t_targ(llm)
142     real :: rh_targ(llm)
143     real :: u_targ(llm)
144     real :: v_targ(llm)
145!
146!---------------------------------------------------------------------
147!  Declarations related to vertical discretization:
148!---------------------------------------------------------------------
149      real :: pzero=1.e5
150      real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
151      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
152
153!---------------------------------------------------------------------
154!  Declarations related to variables
155!---------------------------------------------------------------------
156
157      real :: phi(llm)
158      real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
159      REAL rot(1, llm) ! relative vorticity, in s-1
160      real :: rlat_rad(1),rlon_rad(1)
161      real :: omega(llm),omega2(llm),rho(llm+1)
162      real :: ug(llm),vg(llm),fcoriolis
163      real :: sfdt, cfdt
164      real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
165      real :: w_adv(llm),z_adv(llm)
166      real :: d_t_vert_adv(llm),d_u_vert_adv(llm),d_v_vert_adv(llm)
167      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
168      real :: d_u_nudge(llm),d_v_nudge(llm)
169!      real :: d_u_adv(llm),d_v_adv(llm)
170      real :: d_u_age(llm),d_v_age(llm)
171      real :: alpha
172      real :: ttt
173
174      REAL, ALLOCATABLE, DIMENSION(:,:):: q
175      REAL, ALLOCATABLE, DIMENSION(:,:):: dq
176      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_vert_adv
177      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
178      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
179!      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
180
181!---------------------------------------------------------------------
182!  Initialization of surface variables
183!---------------------------------------------------------------------
184      real :: run_off_lic_0(1)
185      real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
186      real :: tsoil(1,nsoilmx,nbsrf)
187!     real :: agesno(1,nbsrf)
188
189!---------------------------------------------------------------------
190!  Call to phyredem
191!---------------------------------------------------------------------
192      logical :: ok_writedem =.true.
193      real :: sollw_in = 0.
194      real :: solsw_in = 0.
195     
196!---------------------------------------------------------------------
197!  Call to physiq
198!---------------------------------------------------------------------
199      logical :: firstcall=.true.
200      logical :: lastcall=.false.
201      real :: phis(1)    = 0.0
202      real :: dpsrf(1)
203
204!---------------------------------------------------------------------
205!  Initializations of boundary conditions
206!---------------------------------------------------------------------
207      real, allocatable :: phy_nat (:)  ! 0=ocean libre,1=land,2=glacier,3=banquise
208      real, allocatable :: phy_alb (:)  ! Albedo land only (old value condsurf_jyg=0.3)
209      real, allocatable :: phy_sst (:)  ! SST (will not be used; cf read_tsurf1d.F)
210      real, allocatable :: phy_bil (:)  ! Ne sert que pour les slab_ocean
211      real, allocatable :: phy_rug (:) ! Longueur rugosite utilisee sur land only
212      real, allocatable :: phy_ice (:) ! Fraction de glace
213      real, allocatable :: phy_fter(:) ! Fraction de terre
214      real, allocatable :: phy_foce(:) ! Fraction de ocean
215      real, allocatable :: phy_fsic(:) ! Fraction de glace
216      real, allocatable :: phy_flic(:) ! Fraction de glace
217
218!---------------------------------------------------------------------
219!  Fichiers et d'autres variables
220!---------------------------------------------------------------------
221      integer :: k,l,i,it=1,mxcalc
222      integer :: nsrf
223      integer jcode
224      INTEGER read_climoz
225!
226      integer :: it_end ! iteration number of the last call
227!Al1,plev,play,phi,phis,presnivs,
228      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
229      data ecrit_slab_oc/-1/
230!
231!     if flag_inhib_forcing = 0, tendencies of forcing are added
232!                           <> 0, tendencies of forcing are not added
233      INTEGER :: flag_inhib_forcing = 0
234
235
236      print*,'VOUS ENTREZ DANS LE 1D FORMAT STANDARD'
237
238!=====================================================================
239! INITIALIZATIONS
240!=====================================================================
241      du_phys(:)=0.
242      dv_phys(:)=0.
243      dt_phys(:)=0.
244      d_t_vert_adv(:)=0.
245      d_u_vert_adv(:)=0.
246      d_v_vert_adv(:)=0.
247      dt_cooling(:)=0.
248      d_t_adv(:)=0.
249      d_t_nudge(:)=0.
250      d_u_nudge(:)=0.
251      d_v_nudge(:)=0.
252      d_u_adv(:)=0.
253      d_v_adv(:)=0.
254      d_u_age(:)=0.
255      d_v_age(:)=0.
256     
257     
258! Initialization of Common turb_forcing
259       dtime_frcg = 0.
260       Turb_fcg_gcssold=.false.
261       hthturb_gcssold = 0.
262       hqturb_gcssold = 0.
263
264
265
266
267!---------------------------------------------------------------------
268! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
269!---------------------------------------------------------------------
270        call conf_unicol
271!Al1 moves this gcssold var from common fcg_gcssold to
272        Turb_fcg_gcssold = xTurb_fcg_gcssold
273! --------------------------------------------------------------------
274        close(1)
275        write(*,*) 'lmdz1d.def lu => unicol.def'
276
277       forcing_SCM = .true.
278       year_ini_cas=1997
279       ! It is possible that those parameters are run twice.
280       ! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
281
282
283       call getin('anneeref',year_ini_cas)
284       call getin('dayref',day_deb)
285       mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
286       call getin('time_ini',heure_ini_cas)
287
288      print*,'NATURE DE LA SURFACE ',nat_surf
289!
290! Initialization of the logical switch for nudging
291
292     jcode = iflag_nudge
293     do i = 1,nudge_max
294       nudge(i) = mod(jcode,10) .ge. 1
295       jcode = jcode/10
296     enddo
297!-----------------------------------------------------------------------
298!  Definition of the run
299!-----------------------------------------------------------------------
300
301      call conf_gcm( 99, .TRUE. )
302     
303!-----------------------------------------------------------------------
304      allocate( phy_nat (year_len))  ! 0=ocean libre,1=land,2=glacier,3=banquise
305      phy_nat(:)=0.0
306      allocate( phy_alb (year_len))  ! Albedo land only (old value condsurf_jyg=0.3)
307      allocate( phy_sst (year_len))  ! SST (will not be used; cf read_tsurf1d.F)
308      allocate( phy_bil (year_len))  ! Ne sert que pour les slab_ocean
309      phy_bil(:)=1.0
310      allocate( phy_rug (year_len)) ! Longueur rugosite utilisee sur land only
311      allocate( phy_ice (year_len)) ! Fraction de glace
312      phy_ice(:)=0.0
313      allocate( phy_fter(year_len)) ! Fraction de terre
314      phy_fter(:)=0.0
315      allocate( phy_foce(year_len)) ! Fraction de ocean
316      phy_foce(:)=0.0
317      allocate( phy_fsic(year_len)) ! Fraction de glace
318      phy_fsic(:)=0.0
319      allocate( phy_flic(year_len)) ! Fraction de glace
320      phy_flic(:)=0.0
321
322
323!-----------------------------------------------------------------------
324!   Choix du calendrier
325!   -------------------
326
327!      calend = 'earth_365d'
328      if (calend == 'earth_360d') then
329        call ioconf_calendar('360_day')
330        write(*,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
331      else if (calend == 'earth_365d') then
332        call ioconf_calendar('noleap')
333        write(*,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
334      else if (calend == 'earth_366d') then
335        call ioconf_calendar('all_leap')
336        write(*,*)'CALENDRIER CHOISI: Terrestre bissextile'
337      else if (calend == 'gregorian') then
338        stop 'gregorian calend should not be used by normal user'
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
345!-----------------------------------------------------------------------
346!
347!c Date :
348!      La date est supposee donnee sous la forme [annee, numero du jour dans
349!      l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
350!      On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
351!      Le numero du jour est dans "day". L heure est traitee separement.
352!      La date complete est dans "daytime" (l'unite est le jour).
353
354
355      if (nday>0) then
356         fnday=nday
357      else
358         fnday=-nday/float(day_step)
359      endif
360      print *,'fnday=',fnday
361!     start_time doit etre en FRACTION DE JOUR
362      start_time=time_ini/24.
363
364      annee_ref = anneeref
365      mois = 1
366      day_ref = dayref
367      heure = 0.
368      itau_dyn = 0
369      itau_phy = 0
370      call ymds2ju(annee_ref,mois,day_ref,heure,day)
371      day_ini = int(day)
372      day_end = day_ini + int(fnday)
373
374! Convert the initial date to Julian day
375      day_ini_cas=day_deb
376      print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas
377      call ymds2ju                                                         &
378     & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600            &
379     & ,day_ju_ini_cas)
380      print*,'time case 2',day_ini_cas,day_ju_ini_cas
381      daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation
382
383! Print out the actual date of the beginning of the simulation :
384      call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
385      print *,' Time of beginning : ',                                      &
386     &        year_print, month_print, day_print, sec_print
387
388!---------------------------------------------------------------------
389! Initialization of dimensions, geometry and initial state
390!---------------------------------------------------------------------
391!     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
392!     but we still need to initialize dimphy module (klon,klev,etc.)  here.
393      call init_dimphy1D(1,llm)
394      call suphel
395      call init_infotrac
396
397      if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
398      allocate(q(llm,nqtot)) ; q(:,:)=0.
399      allocate(dq(llm,nqtot))
400      allocate(d_q_vert_adv(llm,nqtot))
401      allocate(d_q_adv(llm,nqtot))
402      allocate(d_q_nudge(llm,nqtot))
403!      allocate(d_th_adv(llm))
404
405      q(:,:) = 0.
406      dq(:,:) = 0.
407      d_q_vert_adv(:,:) = 0.
408      d_q_adv(:,:) = 0.
409      d_q_nudge(:,:) = 0.
410
411!
412!   No ozone climatology need be read in this pre-initialization
413!          (phys_state_var_init is called again in physiq)
414      read_climoz = 0
415      nsw=6
416
417      call phys_state_var_init(read_climoz)
418
419      if (ngrid.ne.klon) then
420         print*,'stop in inifis'
421         print*,'Probleme de dimensions :'
422         print*,'ngrid = ',ngrid
423         print*,'klon  = ',klon
424         stop
425      endif
426!!!=====================================================================
427!!! Feedback forcing values for Gateaux differentiation (al1)
428!!!=====================================================================
429!!
430      qsol = qsolinp
431      qsurf = fq_sat(tsurf,psurf/100.)
432      beta_aridity(:,:) = beta_surf
433      day1= day_ini
434      time=daytime-day
435      ts_toga(1)=tsurf ! needed by read_tsurf1d.F
436      rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf))
437
438!
439!! mpl et jyg le 22/08/2012 :
440!!  pour que les cas a flux de surface imposes marchent
441      IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
442       fsens=-wtsurf*rcpd*rho(1)
443       flat=-wqsurf*rlvtt*rho(1)
444       print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
445      ENDIF
446      print*,'Flux sol ',fsens,flat
447
448! Vertical discretization and pressure levels at half and mid levels:
449
450      pa   = 5e4
451!!      preff= 1.01325e5
452      preff = psurf
453      IF (ok_old_disvert) THEN
454        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
455        print *,'On utilise disvert0'
456        aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
457        bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
458        scaleheight=8.
459        pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
460      ELSE
461        call disvert()
462        print *,'On utilise disvert'
463!       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
464!       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
465      ENDIF
466
467      sig_s=presnivs/preff
468      plev =ap+bp*psurf
469      play = 0.5*(plev(1:llm)+plev(2:llm+1))
470      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles.
471
472      IF (forcing_type .eq. 59) THEN
473! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
474      write(*,*) '***********************'
475      do l = 1, llm
476       write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
477       if (trouve_700 .and. play(l).le.70000) then
478         llm700=l
479         print *,'llm700,play=',llm700,play(l)/100.
480         trouve_700= .false.
481       endif
482      enddo
483      write(*,*) '***********************'
484      ENDIF
485
486!
487!=====================================================================
488! EVENTUALLY, READ FORCING DATA :
489!=====================================================================
490
491      INCLUDE "1D_read_forc_cases.h"
492   print*,'A d_t_adv ',d_t_adv(1:20)*86400
493
494      if (forcing_GCM2SCM) then
495        write (*,*) 'forcing_GCM2SCM not yet implemented'
496        stop 'in initialization'
497      endif ! forcing_GCM2SCM
498
499
500!=====================================================================
501! Initialisation de la physique :
502!=====================================================================
503
504!  Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
505!
506! day_step, iphysiq lus dans gcm.def ci-dessus
507! timestep: calcule ci-dessous from rday et day_step
508! ngrid=1
509! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
510! rday: defini dans suphel.F (86400.)
511! day_ini: lu dans run.def (dayref)
512! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
513! airefi,zcufi,zcvfi initialises au debut de ce programme
514! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
515
516
517      day_step = float(nsplit_phys)*day_step/float(iphysiq)
518      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
519      timestep =rday/day_step
520      dtime_frcg = timestep
521!
522      zcufi=airefi
523      zcvfi=airefi
524!
525      rlat_rad(1)=xlat*rpi/180.
526      rlon_rad(1)=xlon*rpi/180.
527
528     ! iniphysiq will call iniaqua who needs year_len from phys_cal_mod
529     year_len_phys_cal_mod=year_len
530           
531     ! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
532     ! e.g. for cell boundaries, which are meaningless in 1D; so pad these
533     ! with '0.' when necessary
534
535      call iniphysiq(iim,jjm,llm, &
536           1,comm_lmdz, &
537           rday,day_ini,timestep,  &
538           (/rlat_rad(1),0./),(/0./), &
539           (/0.,0./),(/rlon_rad(1),0./),  &
540           (/ (/airefi,0./),(/0.,0./) /), &
541           (/zcufi,0.,0.,0./), &
542           (/zcvfi,0./), &
543           ra,rg,rd,rcpd,1)
544      print*,'apres iniphysiq'
545
546! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
547      co2_ppm= 330.0
548      solaire=1370.0
549
550! Ecriture du startphy avant le premier appel a la physique.
551! On le met juste avant pour avoir acces a tous les champs
552
553      if (ok_writedem) then
554
555!--------------------------------------------------------------------------
556! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
557! need : qsol fder snow qsurf evap rugos agesno ftsoil
558!--------------------------------------------------------------------------
559
560        type_ocean = "force"
561        run_off_lic_0(1) = restart_runoff
562        call fonte_neige_init(run_off_lic_0)
563
564        fder=0.
565        snsrf(1,:)=snowmass ! masse de neige des sous surface
566        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
567        fevap=0.
568        z0m(1,:)=rugos     ! couverture de neige des sous surface
569        z0h(1,:)=rugosh    ! couverture de neige des sous surface
570        agesno  = xagesno
571        tsoil(:,:,:)=tsurf
572!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
573!       tsoil(1,1,1)=299.18
574!       tsoil(1,2,1)=300.08
575!       tsoil(1,3,1)=301.88
576!       tsoil(1,4,1)=305.48
577!       tsoil(1,5,1)=308.00
578!       tsoil(1,6,1)=308.00
579!       tsoil(1,7,1)=308.00
580!       tsoil(1,8,1)=308.00
581!       tsoil(1,9,1)=308.00
582!       tsoil(1,10,1)=308.00
583!       tsoil(1,11,1)=308.00
584!-----------------------------------------------------------------------
585        call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
586
587!------------------ prepare limit conditions for limit.nc -----------------
588!--   Ocean force
589
590        print*,'avant phyredem'
591        pctsrf(1,:)=0.
592          if (nat_surf.eq.0.) then
593          pctsrf(1,is_oce)=1.
594          pctsrf(1,is_ter)=0.
595          pctsrf(1,is_lic)=0.
596          pctsrf(1,is_sic)=0.
597        else if (nat_surf .eq. 1) then
598          pctsrf(1,is_oce)=0.
599          pctsrf(1,is_ter)=1.
600          pctsrf(1,is_lic)=0.
601          pctsrf(1,is_sic)=0.
602        else if (nat_surf .eq. 2) then
603          pctsrf(1,is_oce)=0.
604          pctsrf(1,is_ter)=0.
605          pctsrf(1,is_lic)=1.
606          pctsrf(1,is_sic)=0.
607        else if (nat_surf .eq. 3) then
608          pctsrf(1,is_oce)=0.
609          pctsrf(1,is_ter)=0.
610          pctsrf(1,is_lic)=0.
611          pctsrf(1,is_sic)=1.
612
613     end if
614
615
616        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
617     &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
618
619        zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
620        zpic = zpicinp
621        ftsol=tsurf
622        falb_dir=albedo
623        falb_dif=albedo
624        rugoro=rugos
625        t_ancien(1,:)=temp(:)
626        q_ancien(1,:)=q(:,1)
627        ql_ancien = 0.
628        qs_ancien = 0.
629        prlw_ancien = 0.
630        prsw_ancien = 0.
631        prw_ancien = 0.
632!jyg<
633! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
634!!      pbl_tke(:,:,:)=1.e-8
635!        pbl_tke(:,:,:)=0.
636!        pbl_tke(:,2,:)=1.e-2
637!>jyg
638        rain_fall=0.
639        snow_fall=0.
640        solsw=0.
641        solswfdiff=0.
642        sollw=0.
643        sollwdown=rsigma*tsurf**4
644        radsol=0.
645        rnebcon=0.
646        ratqs=0.
647        clwcon=0.
648        zmax0 = 0.
649        zmea=zsurf
650        zstd=0.
651        zsig=0.
652        zgam=0.
653        zval=0.
654        zthe=0.
655        sig1=0.
656        w01=0.
657!
658        wake_deltaq = 0.
659        wake_deltat = 0.
660        wake_delta_pbl_TKE(:,:,:) = 0.
661        delta_tsurf = 0.
662        wake_fip = 0.
663        wake_pe = 0.
664        wake_s = 0.
665        wake_dens = 0.
666        awake_dens = 0.
667        cv_gen = 0.
668        wake_cstar = 0.
669        ale_bl = 0.
670        ale_bl_trig = 0.
671        alp_bl = 0.
672        IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
673        IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
674        entr_therm = 0.
675        detr_therm = 0.
676        f0 = 0.
677        fm_therm = 0.
678        u_ancien(1,:)=u(:)
679        v_ancien(1,:)=v(:)
680        rneb_ancien(1,:)=0.
681 
682        u10m=0.
683        v10m=0.
684        ale_wake=0.
685        ale_bl_stat=0.
686
687!------------------------------------------------------------------------
688! Make file containing restart for the physics (startphy.nc)
689!
690! NB: List of the variables to be written by phyredem (via put_field):
691! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
692! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
693! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
694! radsol,solsw,solswfdiff,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
695! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
696! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
697! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
698! wake_deltat,wake_deltaq,wake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
699! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
700!
701! NB2: The content of the startphy.nc file depends on some flags defined in
702! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
703! to be set at some arbitratry convenient values.
704!------------------------------------------------------------------------
705!Al1 =============== restart option ======================================
706        iflag_physiq=0
707        call getin('iflag_physiq',iflag_physiq)
708
709        if (.not.restart) then
710          iflag_pbl = 5
711          call phyredem ("startphy.nc")
712        else
713! (desallocations)
714        print*,'callin surf final'
715          call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
716        print*,'after surf final'
717          CALL fonte_neige_final(run_off_lic_0)
718        endif
719
720        ok_writedem=.false.
721        print*,'apres phyredem'
722
723      endif ! ok_writedem
724     
725!------------------------------------------------------------------------
726! Make file containing boundary conditions (limit.nc) **Al1->restartdyn***
727! --------------------------------------------------
728! NB: List of the variables to be written in limit.nc
729!     (by writelim.F, subroutine of 1DUTILS.h):
730!        phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
731!        phy_fter,phy_foce,phy_flic,phy_fsic)
732!------------------------------------------------------------------------
733      do i=1,year_len
734        phy_nat(i)  = nat_surf
735        phy_alb(i)  = albedo
736        phy_sst(i)  = tsurf ! read_tsurf1d will be used instead
737        phy_rug(i)  = rugos
738        phy_fter(i) = pctsrf(1,is_ter)
739        phy_foce(i) = pctsrf(1,is_oce)
740        phy_fsic(i) = pctsrf(1,is_sic)
741        phy_flic(i) = pctsrf(1,is_lic)
742      enddo
743
744! fabrication de limit.nc
745      call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,             &
746     &               phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
747
748
749      call phys_state_var_end
750!Al1
751      if (restart) then
752        print*,'call to restart dyn 1d'
753        Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs,          &
754     &              u,v,temp,q,omega2)
755
756       print*,'fnday,annee_ref,day_ref,day_ini',                            &
757     &     fnday,annee_ref,day_ref,day_ini
758!**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
759       day = day_ini
760       day_end = day_ini + nday
761       daytime = day + time_ini/24. ! 1st day and initial time of the simulation
762
763! Print out the actual date of the beginning of the simulation :
764       call ju2ymds(daytime, an, mois, jour, heure)
765       print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
766
767       day = int(daytime)
768       time=daytime-day
769 
770       print*,'****** intialised fields from restart1dyn *******'
771       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
772       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
773       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
774! raz for safety
775       do l=1,llm
776         d_q_vert_adv(l,1) = 0.
777       enddo
778      endif
779!======================  end restart =================================
780      IF (ecrit_slab_oc.eq.1) then
781         open(97,file='div_slab.dat',STATUS='UNKNOWN')
782       elseif (ecrit_slab_oc.eq.0) then
783         open(97,file='div_slab.dat',STATUS='OLD')
784       endif
785!
786!=====================================================================
787#ifdef OUTPUT_PHYS_SCM
788       CALL iophys_ini(timestep)
789#endif
790
791!=====================================================================
792! START OF THE TEMPORAL LOOP :
793!=====================================================================
794           
795      it_end = nint(fnday*day_step)
796      do while(it.le.it_end)
797
798       if (prt_level.ge.1) then
799         print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',                       &
800     &             it,day,time,it_end,day_step
801         print*,'PAS DE TEMPS ',timestep
802       endif
803       if (it.eq.it_end) lastcall=.True.
804
805!---------------------------------------------------------------------
806! Interpolation of forcings in time and onto model levels
807!---------------------------------------------------------------------
808
809      INCLUDE "1D_interp_cases.h"
810
811!---------------------------------------------------------------------
812!  Geopotential :
813!---------------------------------------------------------------------
814        phis(1)=zsurf*RG
815!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
816
817        ! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
818        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
819
820        do l = 1, llm-1
821          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
822     &    (play(l)-play(l+1))/(play(l)+play(l+1))
823        enddo
824
825!---------------------------------------------------------------------
826!  Vertical advection
827!---------------------------------------------------------------------
828
829   IF ( forc_w+forc_omega > 0 ) THEN
830
831      IF ( forc_w == 1 ) THEN
832         w_adv=w_mod_cas
833         z_adv=phi/RG
834      ELSE
835         w_adv=omega
836         z_adv=play
837      ENDIF
838
839      teta=temp*(pzero/play)**rkappa
840      do l=2,llm-1
841        ! vertical tendencies computed as d X / d t = -W  d X / d z
842        d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
843        d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
844        ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
845        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
846        d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))
847      enddo
848      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
849      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
850      d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
851      d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
852   
853   ENDIF
854
855!---------------------------------------------------------------------
856! Listing output for debug prt_level>=1
857!---------------------------------------------------------------------
858       if (prt_level>=1) then
859         print *,' avant physiq : -------- day time ',day,time
860         write(*,*) 'firstcall,lastcall,phis',                               &
861     &               firstcall,lastcall,phis
862       end if
863       if (prt_level>=5) then
864         write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',                   &
865     &        'presniv','plev','play','phi'
866         write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,                   &
867     &         presnivs(l),plev(l),play(l),phi(l),l=1,llm)
868         write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',                    &
869     &         'presniv','u','v','temp','q1','q2','omega2'
870         write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,         &
871     &   presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
872       endif
873
874!---------------------------------------------------------------------
875!   Call physiq :
876!---------------------------------------------------------------------
877       call physiq(ngrid,llm, &
878                    firstcall,lastcall,timestep, &
879                    plev,play,phi,phis,presnivs, &
880                    u,v, rot, temp,q,omega2, &
881                    du_phys,dv_phys,dt_phys,dq,dpsrf)
882                firstcall=.false.
883
884!---------------------------------------------------------------------
885! Listing output for debug
886!---------------------------------------------------------------------
887        if (prt_level>=5) then
888          write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',                  &
889     &        'presniv','plev','play','phi'
890          write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,                  &
891     &    presnivs(l),plev(l),play(l),phi(l),l=1,llm)
892          write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',                   &
893     &         'presniv','u','v','temp','q1','q2','omega2'
894          write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,       &
895     &    presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
896          write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',                   &
897     &         'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'   
898           write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,            &
899     &      presnivs(l),86400*du_phys(l),86400*dv_phys(l),                   &
900     &       86400*dt_phys(l),86400*dq(l,1),dq(l,2),l=1,llm)
901          write(*,*) 'dpsrf',dpsrf
902        endif
903!---------------------------------------------------------------------
904!   Add physical tendencies :
905!---------------------------------------------------------------------
906
907       fcoriolis=2.*sin(rpi*xlat/180.)*romega
908
909      IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
910                                   fcoriolis, xlat,mxcalc
911
912!---------------------------------------------------------------------
913! Geostrophic forcing
914!---------------------------------------------------------------------
915
916      IF ( forc_geo == 0 ) THEN
917              d_u_age(1:mxcalc)=0.
918              d_v_age(1:mxcalc)=0.
919      ELSE
920       sfdt = sin(0.5*fcoriolis*timestep)
921       cfdt = cos(0.5*fcoriolis*timestep)
922
923        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
924     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
925     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
926!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
927!
928       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
929     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
930     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
931!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
932      ENDIF
933!
934!---------------------------------------------------------------------
935!  Nudging
936!  EV: rewrite the section to account for a time- and height-varying
937!  nudging
938!---------------------------------------------------------------------
939      d_t_nudge(:) = 0.
940      d_u_nudge(:) = 0.
941      d_v_nudge(:) = 0.
942      d_q_nudge(:,:) = 0.
943
944      DO l=1,llm
945
946         IF (nudging_u .LT. 0) THEN
947             
948             d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))*invtau_u_nudg_mod_cas(l)
949       
950         ELSE
951
952             IF ( play(l) < p_nudging_u .AND. nint(nudging_u) /= 0 ) &
953             & d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))/nudging_u
954
955         ENDIF
956
957
958         IF (nudging_v .LT. 0) THEN
959             
960             d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))*invtau_v_nudg_mod_cas(l)
961       
962         ELSE
963
964
965             IF ( play(l) < p_nudging_v .AND. nint(nudging_v) /= 0 ) &
966             & d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))/nudging_v
967
968         ENDIF
969
970
971         IF (nudging_t .LT. 0) THEN
972             
973             d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))*invtau_temp_nudg_mod_cas(l)
974       
975         ELSE
976
977
978             IF ( play(l) < p_nudging_t .AND. nint(nudging_t) /= 0 ) &
979             & d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))/nudging_t
980
981          ENDIF
982
983
984         IF (nudging_qv .LT. 0) THEN
985             
986             d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))*invtau_qv_nudg_mod_cas(l)
987       
988         ELSE
989
990             IF ( play(l) < p_nudging_qv .AND. nint(nudging_qv) /= 0 ) &
991             & d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv
992
993         ENDIF
994
995      ENDDO
996
997!---------------------------------------------------------------------
998!  Optional outputs
999!---------------------------------------------------------------------
1000
1001#ifdef OUTPUT_PHYS_SCM
1002      CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
1003      CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv)
1004      CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv)
1005      CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv)
1006      CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1))
1007      CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas)
1008      CALL iophys_ecrit('u',klev,'u','m/s',u)
1009      CALL iophys_ecrit('unud',klev,'unud','m/s',u_nudg_mod_cas)
1010      CALL iophys_ecrit('v',klev,'v','m/s',v)
1011      CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas)
1012      CALL iophys_ecrit('temp',klev,'temp','K',temp)
1013      CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas)
1014      CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge)
1015      CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',1000*86400*d_q_nudge(:,1))
1016#endif
1017
1018    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
1019
1020        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
1021     &              du_phys(1:mxcalc)                                       &
1022     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
1023     &             +d_u_nudge(1:mxcalc) )           
1024        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
1025     &              dv_phys(1:mxcalc)                                       &
1026     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
1027     &             +d_v_nudge(1:mxcalc) )
1028        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
1029     &                dq(1:mxcalc,:)                                        &
1030     &               +d_q_adv(1:mxcalc,:)                                   &
1031     &               +d_q_nudge(1:mxcalc,:) )
1032
1033        if (prt_level.ge.3) then
1034          print *,                                                          &
1035     &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
1036     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
1037           print* ,'dv_phys=',dv_phys
1038           print* ,'d_v_age=',d_v_age
1039           print* ,'d_v_adv=',d_v_adv
1040           print* ,'d_v_nudge=',d_v_nudge
1041           print*, v
1042           print*, vg
1043        endif
1044
1045        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
1046     &              dt_phys(1:mxcalc)                                       &
1047     &             +d_t_adv(1:mxcalc)                                       &
1048     &             +d_t_nudge(1:mxcalc)                                     &
1049     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
1050
1051
1052!=======================================================================
1053!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
1054!=======================================================================
1055
1056        teta=temp*(pzero/play)**rkappa
1057
1058!---------------------------------------------------------------------
1059!   Nudge soil temperature if requested
1060!---------------------------------------------------------------------
1061
1062      IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1063       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
1064     &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
1065      ENDIF
1066
1067!---------------------------------------------------------------------
1068!   Add large-scale tendencies (advection, etc) :
1069!---------------------------------------------------------------------
1070
1071!cc nrlmd
1072!cc        tmpvar=teta
1073!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1074!cc
1075!cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
1076!cc        tmpvar(:)=q(:,1)
1077!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1078!cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
1079!cc        tmpvar(:)=q(:,2)
1080!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1081!cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
1082
1083   END IF ! end if tendency of tendency should be added
1084
1085!---------------------------------------------------------------------
1086!   Air temperature :
1087!---------------------------------------------------------------------       
1088        if (lastcall) then
1089          print*,'Pas de temps final ',it
1090          call ju2ymds(daytime, an, mois, jour, heure)
1091          print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1092        endif
1093
1094!  incremente day time
1095        daytime = daytime+1./day_step
1096        day = int(daytime+0.1/day_step)
1097!        time = max(daytime-day,0.0)
1098!Al1&jyg: correction de bug
1099!cc        time = real(mod(it,day_step))/day_step
1100        time = time_ini/24.+real(mod(it,day_step))/day_step
1101        it=it+1
1102
1103      enddo
1104
1105      if (ecrit_slab_oc.ne.-1) close(97)
1106
1107!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1108! ---------------------------------------------------------------------------
1109       call dyn1dredem("restart1dyn.nc",                                    &
1110     &              plev,play,phi,phis,presnivs,                            &
1111     &              u,v,temp,q,omega2)
1112
1113        CALL abort_gcm ('lmdz1d   ','The End  ',0)
1114
1115END SUBROUTINE scm
Note: See TracBrowser for help on using the repository browser.