source: LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90 @ 4241

Last change on this file since 4241 was 4104, checked in by evignon, 2 years ago

modification des routines pour lecture de la v1 des cas 1D
au format standard + corrections de bug pour le geopotentiel à la surface

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, &
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('360d')
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 infotrac_init
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      print*,'mxcalc=',mxcalc
500!     print*,'zlay=',zlay(mxcalc)
501!      print*,'play=',play(mxcalc)
502
503!! When surface temperature is forced
504      tg= tsurf ! surface T used in read_tsurf1d
505
506
507!=====================================================================
508! Initialisation de la physique :
509!=====================================================================
510
511!  Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
512!
513! day_step, iphysiq lus dans gcm.def ci-dessus
514! timestep: calcule ci-dessous from rday et day_step
515! ngrid=1
516! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
517! rday: defini dans suphel.F (86400.)
518! day_ini: lu dans run.def (dayref)
519! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
520! airefi,zcufi,zcvfi initialises au debut de ce programme
521! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
522
523
524      day_step = float(nsplit_phys)*day_step/float(iphysiq)
525      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
526      timestep =rday/day_step
527      dtime_frcg = timestep
528!
529      zcufi=airefi
530      zcvfi=airefi
531!
532      rlat_rad(1)=xlat*rpi/180.
533      rlon_rad(1)=xlon*rpi/180.
534
535     ! iniphysiq will call iniaqua who needs year_len from phys_cal_mod
536     year_len_phys_cal_mod=year_len
537           
538     ! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
539     ! e.g. for cell boundaries, which are meaningless in 1D; so pad these
540     ! with '0.' when necessary
541
542      call iniphysiq(iim,jjm,llm, &
543           1,comm_lmdz, &
544           rday,day_ini,timestep,  &
545           (/rlat_rad(1),0./),(/0./), &
546           (/0.,0./),(/rlon_rad(1),0./),  &
547           (/ (/airefi,0./),(/0.,0./) /), &
548           (/zcufi,0.,0.,0./), &
549           (/zcvfi,0./), &
550           ra,rg,rd,rcpd,1)
551      print*,'apres iniphysiq'
552
553! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
554      co2_ppm= 330.0
555      solaire=1370.0
556
557! Ecriture du startphy avant le premier appel a la physique.
558! On le met juste avant pour avoir acces a tous les champs
559
560      if (ok_writedem) then
561
562!--------------------------------------------------------------------------
563! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
564! need : qsol fder snow qsurf evap rugos agesno ftsoil
565!--------------------------------------------------------------------------
566
567        type_ocean = "force"
568        run_off_lic_0(1) = restart_runoff
569        call fonte_neige_init(run_off_lic_0)
570
571        fder=0.
572        snsrf(1,:)=snowmass ! masse de neige des sous surface
573        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
574        fevap=0.
575        z0m(1,:)=rugos     ! couverture de neige des sous surface
576        z0h(1,:)=rugosh    ! couverture de neige des sous surface
577        agesno  = xagesno
578        tsoil(:,:,:)=tsurf
579!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
580!       tsoil(1,1,1)=299.18
581!       tsoil(1,2,1)=300.08
582!       tsoil(1,3,1)=301.88
583!       tsoil(1,4,1)=305.48
584!       tsoil(1,5,1)=308.00
585!       tsoil(1,6,1)=308.00
586!       tsoil(1,7,1)=308.00
587!       tsoil(1,8,1)=308.00
588!       tsoil(1,9,1)=308.00
589!       tsoil(1,10,1)=308.00
590!       tsoil(1,11,1)=308.00
591!-----------------------------------------------------------------------
592        call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
593
594!------------------ prepare limit conditions for limit.nc -----------------
595!--   Ocean force
596
597        print*,'avant phyredem'
598        pctsrf(1,:)=0.
599          if (nat_surf.eq.0.) then
600          pctsrf(1,is_oce)=1.
601          pctsrf(1,is_ter)=0.
602          pctsrf(1,is_lic)=0.
603          pctsrf(1,is_sic)=0.
604        else if (nat_surf .eq. 1) then
605          pctsrf(1,is_oce)=0.
606          pctsrf(1,is_ter)=1.
607          pctsrf(1,is_lic)=0.
608          pctsrf(1,is_sic)=0.
609        else if (nat_surf .eq. 2) then
610          pctsrf(1,is_oce)=0.
611          pctsrf(1,is_ter)=0.
612          pctsrf(1,is_lic)=1.
613          pctsrf(1,is_sic)=0.
614        else if (nat_surf .eq. 3) then
615          pctsrf(1,is_oce)=0.
616          pctsrf(1,is_ter)=0.
617          pctsrf(1,is_lic)=0.
618          pctsrf(1,is_sic)=1.
619
620     end if
621
622
623        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
624     &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
625
626        zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
627        zpic = zpicinp
628        ftsol=tsurf
629        falb_dir=albedo
630        falb_dif=albedo
631        rugoro=rugos
632        t_ancien(1,:)=temp(:)
633        q_ancien(1,:)=q(:,1)
634        ql_ancien = 0.
635        qs_ancien = 0.
636        prlw_ancien = 0.
637        prsw_ancien = 0.
638        prw_ancien = 0.
639!jyg<
640! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
641!!      pbl_tke(:,:,:)=1.e-8
642!        pbl_tke(:,:,:)=0.
643!        pbl_tke(:,2,:)=1.e-2
644!>jyg
645        rain_fall=0.
646        snow_fall=0.
647        solsw=0.
648        solswfdiff=0.
649        sollw=0.
650        sollwdown=rsigma*tsurf**4
651        radsol=0.
652        rnebcon=0.
653        ratqs=0.
654        clwcon=0.
655        zmax0 = 0.
656        zmea=zsurf
657        zstd=0.
658        zsig=0.
659        zgam=0.
660        zval=0.
661        zthe=0.
662        sig1=0.
663        w01=0.
664!
665        wake_deltaq = 0.
666        wake_deltat = 0.
667        wake_delta_pbl_TKE(:,:,:) = 0.
668        delta_tsurf = 0.
669        wake_fip = 0.
670        wake_pe = 0.
671        wake_s = 0.
672        wake_dens = 0.
673        awake_dens = 0.
674        cv_gen = 0.
675        wake_cstar = 0.
676        ale_bl = 0.
677        ale_bl_trig = 0.
678        alp_bl = 0.
679        IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
680        IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
681        entr_therm = 0.
682        detr_therm = 0.
683        f0 = 0.
684        fm_therm = 0.
685        u_ancien(1,:)=u(:)
686        v_ancien(1,:)=v(:)
687 
688        u10m=0.
689        v10m=0.
690        ale_wake=0.
691        ale_bl_stat=0.
692
693!------------------------------------------------------------------------
694! Make file containing restart for the physics (startphy.nc)
695!
696! NB: List of the variables to be written by phyredem (via put_field):
697! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
698! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
699! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
700! radsol,solsw,solswfdiff,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
701! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
702! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
703! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
704! wake_deltat,wake_deltaq,wake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
705! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
706!
707! NB2: The content of the startphy.nc file depends on some flags defined in
708! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
709! to be set at some arbitratry convenient values.
710!------------------------------------------------------------------------
711!Al1 =============== restart option ======================================
712        if (.not.restart) then
713          iflag_pbl = 5
714          call phyredem ("startphy.nc")
715        else
716! (desallocations)
717        print*,'callin surf final'
718          call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
719        print*,'after surf final'
720          CALL fonte_neige_final(run_off_lic_0)
721        endif
722
723        ok_writedem=.false.
724        print*,'apres phyredem'
725
726      endif ! ok_writedem
727     
728!------------------------------------------------------------------------
729! Make file containing boundary conditions (limit.nc) **Al1->restartdyn***
730! --------------------------------------------------
731! NB: List of the variables to be written in limit.nc
732!     (by writelim.F, subroutine of 1DUTILS.h):
733!        phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
734!        phy_fter,phy_foce,phy_flic,phy_fsic)
735!------------------------------------------------------------------------
736      do i=1,year_len
737        phy_nat(i)  = nat_surf
738        phy_alb(i)  = albedo
739        phy_sst(i)  = tsurf ! read_tsurf1d will be used instead
740        phy_rug(i)  = rugos
741        phy_fter(i) = pctsrf(1,is_ter)
742        phy_foce(i) = pctsrf(1,is_oce)
743        phy_fsic(i) = pctsrf(1,is_sic)
744        phy_flic(i) = pctsrf(1,is_lic)
745      enddo
746
747! fabrication de limit.nc
748      call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,             &
749     &               phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
750
751
752      call phys_state_var_end
753!Al1
754      if (restart) then
755        print*,'call to restart dyn 1d'
756        Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs,          &
757     &              u,v,temp,q,omega2)
758
759       print*,'fnday,annee_ref,day_ref,day_ini',                            &
760     &     fnday,annee_ref,day_ref,day_ini
761!**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
762       day = day_ini
763       day_end = day_ini + nday
764       daytime = day + time_ini/24. ! 1st day and initial time of the simulation
765
766! Print out the actual date of the beginning of the simulation :
767       call ju2ymds(daytime, an, mois, jour, heure)
768       print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
769
770       day = int(daytime)
771       time=daytime-day
772 
773       print*,'****** intialised fields from restart1dyn *******'
774       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
775       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
776       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
777! raz for safety
778       do l=1,llm
779         d_q_vert_adv(l,1) = 0.
780       enddo
781      endif
782!======================  end restart =================================
783      IF (ecrit_slab_oc.eq.1) then
784         open(97,file='div_slab.dat',STATUS='UNKNOWN')
785       elseif (ecrit_slab_oc.eq.0) then
786         open(97,file='div_slab.dat',STATUS='OLD')
787       endif
788!
789!=====================================================================
790#ifdef OUTPUT_PHYS_SCM
791       CALL iophys_ini(timestep)
792#endif
793
794!=====================================================================
795! START OF THE TEMPORAL LOOP :
796!=====================================================================
797           
798      it_end = nint(fnday*day_step)
799      do while(it.le.it_end)
800
801       if (prt_level.ge.1) then
802         print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',                       &
803     &             it,day,time,it_end,day_step
804         print*,'PAS DE TEMPS ',timestep
805       endif
806       if (it.eq.it_end) lastcall=.True.
807
808!---------------------------------------------------------------------
809! Interpolation of forcings in time and onto model levels
810!---------------------------------------------------------------------
811
812#include "1D_interp_cases.h"
813
814!---------------------------------------------------------------------
815!  Geopotential :
816!---------------------------------------------------------------------
817        phis(1)=zsurf*RG
818!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
819
820        ! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
821        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
822
823        do l = 1, llm-1
824          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
825     &    (play(l)-play(l+1))/(play(l)+play(l+1))
826        enddo
827
828!---------------------------------------------------------------------
829!  Vertical advection
830!---------------------------------------------------------------------
831
832   IF ( forc_w+forc_omega > 0 ) THEN
833
834      IF ( forc_w == 1 ) THEN
835         w_adv=w_mod_cas
836         z_adv=phi/RG
837      ELSE
838         w_adv=omega
839         z_adv=play
840      ENDIF
841
842      teta=temp*(pzero/play)**rkappa
843      do l=2,llm-1
844        ! vertical tendencies computed as d X / d t = -W  d X / d z
845        d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
846        d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
847        ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
848        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
849        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))
850      enddo
851      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
852      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
853      d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
854      d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
855   
856   ENDIF
857
858!---------------------------------------------------------------------
859! Listing output for debug prt_level>=1
860!---------------------------------------------------------------------
861       if (prt_level>=1) then
862         print *,' avant physiq : -------- day time ',day,time
863         write(*,*) 'firstcall,lastcall,phis',                               &
864     &               firstcall,lastcall,phis
865       end if
866       if (prt_level>=5) then
867         write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',                   &
868     &        'presniv','plev','play','phi'
869         write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,                   &
870     &         presnivs(l),plev(l),play(l),phi(l),l=1,llm)
871         write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',                    &
872     &         'presniv','u','v','temp','q1','q2','omega2'
873         write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,         &
874     &   presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
875       endif
876
877!---------------------------------------------------------------------
878!   Call physiq :
879!---------------------------------------------------------------------
880       call physiq(ngrid,llm, &
881                    firstcall,lastcall,timestep, &
882                    plev,play,phi,phis,presnivs, &
883                    u,v, rot, temp,q,omega2, &
884                    du_phys,dv_phys,dt_phys,dq,dpsrf)
885                firstcall=.false.
886
887!---------------------------------------------------------------------
888! Listing output for debug
889!---------------------------------------------------------------------
890        if (prt_level>=5) then
891          write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',                  &
892     &        'presniv','plev','play','phi'
893          write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,                  &
894     &    presnivs(l),plev(l),play(l),phi(l),l=1,llm)
895          write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',                   &
896     &         'presniv','u','v','temp','q1','q2','omega2'
897          write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,       &
898     &    presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
899          write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',                   &
900     &         'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'   
901           write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,            &
902     &      presnivs(l),86400*du_phys(l),86400*dv_phys(l),                   &
903     &       86400*dt_phys(l),86400*dq(l,1),dq(l,2),l=1,llm)
904          write(*,*) 'dpsrf',dpsrf
905        endif
906!---------------------------------------------------------------------
907!   Add physical tendencies :
908!---------------------------------------------------------------------
909
910       fcoriolis=2.*sin(rpi*xlat/180.)*romega
911
912      IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
913                                   fcoriolis, xlat,mxcalc
914
915!---------------------------------------------------------------------
916! Geostrophic forcing
917!---------------------------------------------------------------------
918
919      IF ( forc_geo == 0 ) THEN
920              d_u_age(1:mxcalc)=0.
921              d_v_age(1:mxcalc)=0.
922      ELSE
923       sfdt = sin(0.5*fcoriolis*timestep)
924       cfdt = cos(0.5*fcoriolis*timestep)
925
926        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
927     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
928     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
929!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
930!
931       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
932     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
933     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
934!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
935      ENDIF
936!
937!---------------------------------------------------------------------
938!  Nudging
939!  EV: rewrite the section to account for a time- and height-varying
940!  nudging
941!---------------------------------------------------------------------
942      d_t_nudge(:) = 0.
943      d_u_nudge(:) = 0.
944      d_v_nudge(:) = 0.
945      d_q_nudge(:,:) = 0.
946
947      DO l=1,llm
948
949         IF (nudging_u .LT. 0) THEN
950             
951             d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))*invtau_u_nudg_mod_cas(l)
952       
953         ELSE
954
955             IF ( play(l) < p_nudging_u .AND. nint(nudging_u) /= 0 ) &
956             & d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))/nudging_u
957
958         ENDIF
959
960
961         IF (nudging_v .LT. 0) THEN
962             
963             d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))*invtau_v_nudg_mod_cas(l)
964       
965         ELSE
966
967
968             IF ( play(l) < p_nudging_v .AND. nint(nudging_v) /= 0 ) &
969             & d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))/nudging_v
970
971         ENDIF
972
973
974         IF (nudging_t .LT. 0) THEN
975             
976             d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))*invtau_temp_nudg_mod_cas(l)
977       
978         ELSE
979
980
981             IF ( play(l) < p_nudging_t .AND. nint(nudging_t) /= 0 ) &
982             & d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))/nudging_t
983
984          ENDIF
985
986
987         IF (nudging_qv .LT. 0) THEN
988             
989             d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))*invtau_qv_nudg_mod_cas(l)
990       
991         ELSE
992
993             IF ( play(l) < p_nudging_qv .AND. nint(nudging_qv) /= 0 ) &
994             & d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv
995
996         ENDIF
997
998      ENDDO
999
1000!---------------------------------------------------------------------
1001!  Optional outputs
1002!---------------------------------------------------------------------
1003
1004#ifdef OUTPUT_PHYS_SCM
1005      CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
1006      CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv)
1007      CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv)
1008      CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv)
1009      CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1))
1010      CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas)
1011      CALL iophys_ecrit('u',klev,'u','m/s',u)
1012      CALL iophys_ecrit('unud',klev,'unud','m/s',u_nudg_mod_cas)
1013      CALL iophys_ecrit('v',klev,'v','m/s',v)
1014      CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas)
1015      CALL iophys_ecrit('temp',klev,'temp','K',temp)
1016      CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas)
1017      CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge)
1018      CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',1000*86400*d_q_nudge(:,1))
1019#endif
1020
1021    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
1022
1023        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
1024     &              du_phys(1:mxcalc)                                       &
1025     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
1026     &             +d_u_nudge(1:mxcalc) )           
1027        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
1028     &              dv_phys(1:mxcalc)                                       &
1029     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
1030     &             +d_v_nudge(1:mxcalc) )
1031        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
1032     &                dq(1:mxcalc,:)                                        &
1033     &               +d_q_adv(1:mxcalc,:)                                   &
1034     &               +d_q_nudge(1:mxcalc,:) )
1035
1036        if (prt_level.ge.3) then
1037          print *,                                                          &
1038     &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
1039     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
1040           print* ,'dv_phys=',dv_phys
1041           print* ,'d_v_age=',d_v_age
1042           print* ,'d_v_adv=',d_v_adv
1043           print* ,'d_v_nudge=',d_v_nudge
1044           print*, v
1045           print*, vg
1046        endif
1047
1048        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
1049     &              dt_phys(1:mxcalc)                                       &
1050     &             +d_t_adv(1:mxcalc)                                       &
1051     &             +d_t_nudge(1:mxcalc)                                     &
1052     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
1053
1054
1055!=======================================================================
1056!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
1057!=======================================================================
1058
1059        teta=temp*(pzero/play)**rkappa
1060
1061!---------------------------------------------------------------------
1062!   Nudge soil temperature if requested
1063!---------------------------------------------------------------------
1064
1065      IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1066       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
1067     &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
1068      ENDIF
1069
1070!---------------------------------------------------------------------
1071!   Add large-scale tendencies (advection, etc) :
1072!---------------------------------------------------------------------
1073
1074!cc nrlmd
1075!cc        tmpvar=teta
1076!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1077!cc
1078!cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
1079!cc        tmpvar(:)=q(:,1)
1080!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1081!cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
1082!cc        tmpvar(:)=q(:,2)
1083!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
1084!cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
1085
1086   END IF ! end if tendency of tendency should be added
1087
1088!---------------------------------------------------------------------
1089!   Air temperature :
1090!---------------------------------------------------------------------       
1091        if (lastcall) then
1092          print*,'Pas de temps final ',it
1093          call ju2ymds(daytime, an, mois, jour, heure)
1094          print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1095        endif
1096
1097!  incremente day time
1098        daytime = daytime+1./day_step
1099        day = int(daytime+0.1/day_step)
1100!        time = max(daytime-day,0.0)
1101!Al1&jyg: correction de bug
1102!cc        time = real(mod(it,day_step))/day_step
1103        time = time_ini/24.+real(mod(it,day_step))/day_step
1104        it=it+1
1105
1106      enddo
1107
1108      if (ecrit_slab_oc.ne.-1) close(97)
1109
1110!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1111! ---------------------------------------------------------------------------
1112       call dyn1dredem("restart1dyn.nc",                                    &
1113     &              plev,play,phi,phis,presnivs,                            &
1114     &              u,v,temp,q,omega2)
1115
1116        CALL abort_gcm ('lmdz1d   ','The End  ',0)
1117
1118END SUBROUTINE scm
Note: See TracBrowser for help on using the repository browser.