source: LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 4128

Last change on this file since 4128 was 2990, checked in by fhourdin, 7 years ago

Modif 1D, lecture format standard. Interpolation en temperature potentielle.

File size: 174.3 KB
Line 
1#include "conf_gcm.F90"
2
3!
4! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
5!
6!
7!
8      SUBROUTINE conf_unicol
9!
10#ifdef CPP_IOIPSL
11      use IOIPSL
12#else
13! if not using IOIPSL, we still need to use (a local version of) getin
14      use ioipsl_getincom
15#endif
16      USE print_control_mod, ONLY: lunout
17      IMPLICIT NONE
18!-----------------------------------------------------------------------
19!     Auteurs :   A. Lahellec  .
20!
21!   Declarations :
22!   --------------
23
24#include "compar1d.h"
25#include "flux_arp.h"
26#include "tsoilnudge.h"
27#include "fcg_gcssold.h"
28#include "fcg_racmo.h"
29!
30!
31!   local:
32!   ------
33
34!      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
35     
36!
37!  -------------------------------------------------------------------
38!
39!      .........    Initilisation parametres du lmdz1D      ..........
40!
41!---------------------------------------------------------------------
42!   initialisations:
43!   ----------------
44
45!Config  Key  = lunout
46!Config  Desc = unite de fichier pour les impressions
47!Config  Def  = 6
48!Config  Help = unite de fichier pour les impressions
49!Config         (defaut sortie standard = 6)
50      lunout=6
51!      CALL getin('lunout', lunout)
52      IF (lunout /= 5 .and. lunout /= 6) THEN
53        OPEN(lunout,FILE='lmdz.out')
54      ENDIF
55
56!Config  Key  = prt_level
57!Config  Desc = niveau d'impressions de debogage
58!Config  Def  = 0
59!Config  Help = Niveau d'impression pour le debogage
60!Config         (0 = minimum d'impression)
61!      prt_level = 0
62!      CALL getin('prt_level',prt_level)
63
64!-----------------------------------------------------------------------
65!  Parametres de controle du run:
66!-----------------------------------------------------------------------
67
68!Config  Key  = restart
69!Config  Desc = on repart des startphy et start1dyn
70!Config  Def  = false
71!Config  Help = les fichiers restart doivent etre renomme en start
72       restart =.false.
73       CALL getin('restart',restart)
74
75!Config  Key  = forcing_type
76!Config  Desc = defines the way the SCM is forced:
77!Config  Def  = 0
78!!Config  Help = 0 ==> forcing_les = .true.
79!             initial profiles from file prof.inp.001
80!             no forcing by LS convergence ; 
81!             surface temperature imposed ;
82!             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
83!         = 1 ==> forcing_radconv = .true.
84!             idem forcing_type = 0, but the imposed radiative cooling
85!             is set to 0 (hence, if iflag_radia=0 in physiq.def, 
86!             then there is no radiative cooling at all)
87!         = 2 ==> forcing_toga = .true.
88!             initial profiles from TOGA-COARE IFA files
89!             LS convergence and SST imposed from TOGA-COARE IFA files
90!         = 3 ==> forcing_GCM2SCM = .true.
91!             initial profiles from the GCM output
92!             LS convergence imposed from the GCM output
93!         = 4 ==> forcing_twpi = .true.
94!             initial profiles from TWPICE nc files
95!             LS convergence and SST imposed from TWPICE nc files
96!         = 5 ==> forcing_rico = .true.
97!             initial profiles from RICO idealized
98!             LS convergence imposed from  RICO (cst) 
99!         = 6 ==> forcing_amma = .true.
100!         = 10 ==> forcing_case = .true.
101!             initial profiles from case.nc file
102!         = 40 ==> forcing_GCSSold = .true.
103!             initial profile from GCSS file
104!             LS convergence imposed from GCSS file
105!         = 50 ==> forcing_fire = .true.
106!         = 59 ==> forcing_sandu = .true.
107!             initial profiles from sanduref file: see prof.inp.001
108!             SST varying with time and divergence constante: see ifa_sanduref.txt file
109!             Radiation has to be computed interactively
110!         = 60 ==> forcing_astex = .true.
111!             initial profiles from file: see prof.inp.001 
112!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
113!             Radiation has to be computed interactively
114!         = 61 ==> forcing_armcu = .true.
115!             initial profiles from file: see prof.inp.001 
116!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
117!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
118!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
119!             Radiation to be switched off
120!         > 100 ==> forcing_case = .true. or forcing_case2 = .true.
121!             initial profiles from case.nc file
122!
123       forcing_type = 0
124       CALL getin('forcing_type',forcing_type)
125         imp_fcg_gcssold   = .false.
126         ts_fcg_gcssold    = .false. 
127         Tp_fcg_gcssold    = .false. 
128         Tp_ini_gcssold    = .false. 
129         xTurb_fcg_gcssold = .false. 
130        IF (forcing_type .eq.40) THEN
131          CALL getin('imp_fcg',imp_fcg_gcssold)
132          CALL getin('ts_fcg',ts_fcg_gcssold)
133          CALL getin('tp_fcg',Tp_fcg_gcssold)
134          CALL getin('tp_ini',Tp_ini_gcssold)
135          CALL getin('turb_fcg',xTurb_fcg_gcssold)
136        ENDIF
137
138!Parametres de forcage
139!Config  Key  = tend_t
140!Config  Desc = forcage ou non par advection de T
141!Config  Def  = false
142!Config  Help = forcage ou non par advection de T
143       tend_t =0
144       CALL getin('tend_t',tend_t)
145
146!Config  Key  = tend_q
147!Config  Desc = forcage ou non par advection de q
148!Config  Def  = false
149!Config  Help = forcage ou non par advection de q
150       tend_q =0
151       CALL getin('tend_q',tend_q)
152
153!Config  Key  = tend_u
154!Config  Desc = forcage ou non par advection de u
155!Config  Def  = false
156!Config  Help = forcage ou non par advection de u
157       tend_u =0
158       CALL getin('tend_u',tend_u)
159
160!Config  Key  = tend_v
161!Config  Desc = forcage ou non par advection de v
162!Config  Def  = false
163!Config  Help = forcage ou non par advection de v
164       tend_v =0
165       CALL getin('tend_v',tend_v)
166
167!Config  Key  = tend_w
168!Config  Desc = forcage ou non par vitesse verticale
169!Config  Def  = false
170!Config  Help = forcage ou non par vitesse verticale
171       tend_w =0
172       CALL getin('tend_w',tend_w)
173
174!Config  Key  = tend_rayo
175!Config  Desc = forcage ou non par dtrad
176!Config  Def  = false
177!Config  Help = forcage ou non par dtrad
178       tend_rayo =0
179       CALL getin('tend_rayo',tend_rayo)
180
181
182!Config  Key  = nudge_t
183!Config  Desc = constante de nudging de T
184!Config  Def  = false
185!Config  Help = constante de nudging de T
186       nudge_t =0.
187       CALL getin('nudge_t',nudge_t)
188
189!Config  Key  = nudge_q
190!Config  Desc = constante de nudging de q
191!Config  Def  = false
192!Config  Help = constante de nudging de q
193       nudge_q =0.
194       CALL getin('nudge_q',nudge_q)
195
196!Config  Key  = nudge_u
197!Config  Desc = constante de nudging de u
198!Config  Def  = false
199!Config  Help = constante de nudging de u
200       nudge_u =0.
201       CALL getin('nudge_u',nudge_u)
202
203!Config  Key  = nudge_v
204!Config  Desc = constante de nudging de v
205!Config  Def  = false
206!Config  Help = constante de nudging de v
207       nudge_v =0.
208       CALL getin('nudge_v',nudge_v)
209
210!Config  Key  = nudge_w
211!Config  Desc = constante de nudging de w
212!Config  Def  = false
213!Config  Help = constante de nudging de w
214       nudge_w =0.
215       CALL getin('nudge_w',nudge_w)
216
217
218!Config  Key  = iflag_nudge
219!Config  Desc = atmospheric nudging ttype (decimal code)
220!Config  Def  = 0
221!Config  Help = 0 ==> no nudging
222!  If digit number n of iflag_nudge is set, then nudging of type n is on
223!  If digit number n of iflag_nudge is not set, then nudging of type n is off
224!   (digits are numbered from the right)
225       iflag_nudge = 0
226       CALL getin('iflag_nudge',iflag_nudge)
227
228!Config  Key  = ok_flux_surf
229!Config  Desc = forcage ou non par les flux de surface
230!Config  Def  = false
231!Config  Help = forcage ou non par les flux de surface
232       ok_flux_surf =.false.
233       CALL getin('ok_flux_surf',ok_flux_surf)
234
235!Config  Key  = ok_prescr_ust
236!Config  Desc = ustar impose ou non
237!Config  Def  = false
238!Config  Help = ustar impose ou non
239       ok_prescr_ust = .false.
240       CALL getin('ok_prescr_ust',ok_prescr_ust)
241
242!Config  Key  = ok_old_disvert
243!Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
244!Config  Def  = false
245!Config  Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
246       ok_old_disvert = .FALSE.
247       CALL getin('ok_old_disvert',ok_old_disvert)
248
249!Config  Key  = time_ini
250!Config  Desc = meaningless in this  case
251!Config  Def  = 0.
252!Config  Help = 
253       tsurf = 0.
254       CALL getin('time_ini',time_ini)
255
256!Config  Key  = rlat et rlon
257!Config  Desc = latitude et longitude
258!Config  Def  = 0.0  0.0
259!Config  Help = fixe la position de la colonne
260       xlat = 0.
261       xlon = 0.
262       CALL getin('rlat',xlat)
263       CALL getin('rlon',xlon)
264
265!Config  Key  = airephy
266!Config  Desc = Grid cell area
267!Config  Def  = 1.e11
268!Config  Help = 
269       airefi = 1.e11
270       CALL getin('airephy',airefi)
271
272!Config  Key  = nat_surf
273!Config  Desc = surface type
274!Config  Def  = 0 (ocean)
275!Config  Help = 0=ocean,1=land,2=glacier,3=banquise
276       nat_surf = 0.
277       CALL getin('nat_surf',nat_surf)
278
279!Config  Key  = tsurf
280!Config  Desc = surface temperature
281!Config  Def  = 290.
282!Config  Help = not used if type_ts_forcing=1 in lmdz1d.F
283       tsurf = 290.
284       CALL getin('tsurf',tsurf)
285
286!Config  Key  = psurf
287!Config  Desc = surface pressure
288!Config  Def  = 102400.
289!Config  Help = 
290       psurf = 102400.
291       CALL getin('psurf',psurf)
292
293!Config  Key  = zsurf
294!Config  Desc = surface altitude
295!Config  Def  = 0.
296!Config  Help = 
297       zsurf = 0.
298       CALL getin('zsurf',zsurf)
299
300!Config  Key  = rugos
301!Config  Desc = coefficient de frottement
302!Config  Def  = 0.0001
303!Config  Help = calcul du Cdrag
304       rugos = 0.0001
305       CALL getin('rugos',rugos)
306
307!Config  Key  = rugosh
308!Config  Desc = coefficient de frottement
309!Config  Def  = rugos
310!Config  Help = calcul du Cdrag
311       rugosh = rugos
312       CALL getin('rugosh',rugosh)
313
314
315
316!Config  Key  = snowmass
317!Config  Desc = mass de neige de la surface en kg/m2
318!Config  Def  = 0.0000
319!Config  Help = snowmass
320       snowmass = 0.0000
321       CALL getin('snowmass',snowmass)
322
323!Config  Key  = wtsurf et wqsurf
324!Config  Desc = ???
325!Config  Def  = 0.0 0.0
326!Config  Help = 
327       wtsurf = 0.0
328       wqsurf = 0.0
329       CALL getin('wtsurf',wtsurf)
330       CALL getin('wqsurf',wqsurf)
331
332!Config  Key  = albedo
333!Config  Desc = albedo
334!Config  Def  = 0.09
335!Config  Help = 
336       albedo = 0.09
337       CALL getin('albedo',albedo)
338
339!Config  Key  = agesno
340!Config  Desc = age de la neige
341!Config  Def  = 30.0
342!Config  Help = 
343       xagesno = 30.0
344       CALL getin('agesno',xagesno)
345
346!Config  Key  = restart_runoff
347!Config  Desc = age de la neige
348!Config  Def  = 30.0
349!Config  Help = 
350       restart_runoff = 0.0
351       CALL getin('restart_runoff',restart_runoff)
352
353!Config  Key  = qsolinp
354!Config  Desc = initial bucket water content (kg/m2) when land (5std)
355!Config  Def  = 30.0
356!Config  Help = 
357       qsolinp = 1.
358       CALL getin('qsolinp',qsolinp)
359
360!Config  Key  = zpicinp
361!Config  Desc = denivellation orographie
362!Config  Def  = 0.
363!Config  Help =  input brise
364       zpicinp = 0.
365       CALL getin('zpicinp',zpicinp)
366!Config key = nudge_tsoil
367!Config  Desc = activation of soil temperature nudging
368!Config  Def  = .FALSE.
369!Config  Help = ...
370
371       nudge_tsoil=.FALSE.
372       CALL getin('nudge_tsoil',nudge_tsoil)
373
374!Config key = isoil_nudge
375!Config  Desc = level number where soil temperature is nudged
376!Config  Def  = 3
377!Config  Help = ...
378
379       isoil_nudge=3
380       CALL getin('isoil_nudge',isoil_nudge)
381
382!Config key = Tsoil_nudge
383!Config  Desc = target temperature for tsoil(isoil_nudge)
384!Config  Def  = 300.
385!Config  Help = ...
386
387       Tsoil_nudge=300.
388       CALL getin('Tsoil_nudge',Tsoil_nudge)
389
390!Config key = tau_soil_nudge
391!Config  Desc = nudging relaxation time for tsoil
392!Config  Def  = 3600.
393!Config  Help = ...
394
395       tau_soil_nudge=3600.
396       CALL getin('tau_soil_nudge',tau_soil_nudge)
397
398!----------------------------------------------------------
399! Param??tres de for??age pour les forcages communs:
400! Pour les forcages communs: ces entiers valent 0 ou 1
401! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
402! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale
403! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv)
404! forcages en omega, w, vent geostrophique ou ustar
405! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging
406!----------------------------------------------------------
407
408!Config  Key  = tadv
409!Config  Desc = forcage ou non par advection totale de T
410!Config  Def  = false
411!Config  Help = forcage ou non par advection totale de T
412       tadv =0
413       CALL getin('tadv',tadv)
414
415!Config  Key  = tadvv
416!Config  Desc = forcage ou non par advection verticale de T
417!Config  Def  = false
418!Config  Help = forcage ou non par advection verticale de T
419       tadvv =0
420       CALL getin('tadvv',tadvv)
421
422!Config  Key  = tadvh
423!Config  Desc = forcage ou non par advection horizontale de T
424!Config  Def  = false
425!Config  Help = forcage ou non par advection horizontale de T
426       tadvh =0
427       CALL getin('tadvh',tadvh)
428
429!Config  Key  = thadv
430!Config  Desc = forcage ou non par advection totale de Theta
431!Config  Def  = false
432!Config  Help = forcage ou non par advection totale de Theta
433       thadv =0
434       CALL getin('thadv',thadv)
435
436!Config  Key  = thadvv
437!Config  Desc = forcage ou non par advection verticale de Theta
438!Config  Def  = false
439!Config  Help = forcage ou non par advection verticale de Theta
440       thadvv =0
441       CALL getin('thadvv',thadvv)
442
443!Config  Key  = thadvh
444!Config  Desc = forcage ou non par advection horizontale de Theta
445!Config  Def  = false
446!Config  Help = forcage ou non par advection horizontale de Theta
447       thadvh =0
448       CALL getin('thadvh',thadvh)
449
450!Config  Key  = qadv
451!Config  Desc = forcage ou non par advection totale de Q
452!Config  Def  = false
453!Config  Help = forcage ou non par advection totale de Q
454       qadv =0
455       CALL getin('qadv',qadv)
456
457!Config  Key  = qadvv
458!Config  Desc = forcage ou non par advection verticale de Q
459!Config  Def  = false
460!Config  Help = forcage ou non par advection verticale de Q
461       qadvv =0
462       CALL getin('qadvv',qadvv)
463
464!Config  Key  = qadvh
465!Config  Desc = forcage ou non par advection horizontale de Q
466!Config  Def  = false
467!Config  Help = forcage ou non par advection horizontale de Q
468       qadvh =0
469       CALL getin('qadvh',qadvh)
470
471!Config  Key  = trad
472!Config  Desc = forcage ou non par tendance radiative
473!Config  Def  = false
474!Config  Help = forcage ou non par tendance radiative
475       trad =0
476       CALL getin('trad',trad)
477
478!Config  Key  = forc_omega
479!Config  Desc = forcage ou non par omega
480!Config  Def  = false
481!Config  Help = forcage ou non par omega
482       forc_omega =0
483       CALL getin('forc_omega',forc_omega)
484
485!Config  Key  = forc_u
486!Config  Desc = forcage ou non par u
487!Config  Def  = false
488!Config  Help = forcage ou non par u
489       forc_u =0
490       CALL getin('forc_u',forc_u)
491
492!Config  Key  = forc_v
493!Config  Desc = forcage ou non par v
494!Config  Def  = false
495!Config  Help = forcage ou non par v
496       forc_v =0
497       CALL getin('forc_v',forc_v)
498!Config  Key  = forc_w
499!Config  Desc = forcage ou non par w
500!Config  Def  = false
501!Config  Help = forcage ou non par w
502       forc_w =0
503       CALL getin('forc_w',forc_w)
504
505!Config  Key  = forc_geo
506!Config  Desc = forcage ou non par geo
507!Config  Def  = false
508!Config  Help = forcage ou non par geo
509       forc_geo =0
510       CALL getin('forc_geo',forc_geo)
511
512! Meme chose que ok_precr_ust
513!Config  Key  = forc_ustar
514!Config  Desc = forcage ou non par ustar
515!Config  Def  = false
516!Config  Help = forcage ou non par ustar
517       forc_ustar =0
518       CALL getin('forc_ustar',forc_ustar)
519       IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
520
521!Config  Key  = nudging_u
522!Config  Desc = forcage ou non par nudging sur u
523!Config  Def  = false
524!Config  Help = forcage ou non par nudging sur u
525       nudging_u =0
526       CALL getin('nudging_u',nudging_u)
527
528!Config  Key  = nudging_v
529!Config  Desc = forcage ou non par nudging sur v
530!Config  Def  = false
531!Config  Help = forcage ou non par nudging sur v
532       nudging_v =0
533       CALL getin('nudging_v',nudging_v)
534
535!Config  Key  = nudging_w
536!Config  Desc = forcage ou non par nudging sur w
537!Config  Def  = false
538!Config  Help = forcage ou non par nudging sur w
539       nudging_w =0
540       CALL getin('nudging_w',nudging_w)
541
542!Config  Key  = nudging_q
543!Config  Desc = forcage ou non par nudging sur q
544!Config  Def  = false
545!Config  Help = forcage ou non par nudging sur q
546       nudging_q =0
547       CALL getin('nudging_q',nudging_q)
548
549!Config  Key  = nudging_t
550!Config  Desc = forcage ou non par nudging sur t
551!Config  Def  = false
552!Config  Help = forcage ou non par nudging sur t
553       nudging_t =0
554       CALL getin('nudging_t',nudging_t)
555
556
557
558      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
559      write(lunout,*)' Configuration des parametres du gcm1D: '
560      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
561      write(lunout,*)' restart = ', restart
562      write(lunout,*)' forcing_type = ', forcing_type
563      write(lunout,*)' time_ini = ', time_ini
564      write(lunout,*)' rlat = ', xlat
565      write(lunout,*)' rlon = ', xlon
566      write(lunout,*)' airephy = ', airefi
567      write(lunout,*)' nat_surf = ', nat_surf
568      write(lunout,*)' tsurf = ', tsurf
569      write(lunout,*)' psurf = ', psurf
570      write(lunout,*)' zsurf = ', zsurf
571      write(lunout,*)' rugos = ', rugos
572      write(lunout,*)' snowmass=', snowmass
573      write(lunout,*)' wtsurf = ', wtsurf
574      write(lunout,*)' wqsurf = ', wqsurf
575      write(lunout,*)' albedo = ', albedo
576      write(lunout,*)' xagesno = ', xagesno
577      write(lunout,*)' restart_runoff = ', restart_runoff
578      write(lunout,*)' qsolinp = ', qsolinp
579      write(lunout,*)' zpicinp = ', zpicinp
580      write(lunout,*)' nudge_tsoil = ', nudge_tsoil
581      write(lunout,*)' isoil_nudge = ', isoil_nudge
582      write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
583      write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
584      write(lunout,*)' tadv =      ', tadv
585      write(lunout,*)' tadvv =     ', tadvv
586      write(lunout,*)' tadvh =     ', tadvh
587      write(lunout,*)' thadv =     ', thadv
588      write(lunout,*)' thadvv =    ', thadvv
589      write(lunout,*)' thadvh =    ', thadvh
590      write(lunout,*)' qadv =      ', qadv
591      write(lunout,*)' qadvv =     ', qadvv
592      write(lunout,*)' qadvh =     ', qadvh
593      write(lunout,*)' trad =      ', trad
594      write(lunout,*)' forc_omega = ', forc_omega
595      write(lunout,*)' forc_w     = ', forc_w
596      write(lunout,*)' forc_geo   = ', forc_geo
597      write(lunout,*)' forc_ustar = ', forc_ustar
598      write(lunout,*)' nudging_u  = ', nudging_u
599      write(lunout,*)' nudging_v  = ', nudging_v
600      write(lunout,*)' nudging_t  = ', nudging_t
601      write(lunout,*)' nudging_q  = ', nudging_q
602      IF (forcing_type .eq.40) THEN
603        write(lunout,*) '--- Forcing type GCSS Old --- with:'
604        write(lunout,*)'imp_fcg',imp_fcg_gcssold
605        write(lunout,*)'ts_fcg',ts_fcg_gcssold
606        write(lunout,*)'tp_fcg',Tp_fcg_gcssold
607        write(lunout,*)'tp_ini',Tp_ini_gcssold
608        write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
609      ENDIF
610
611      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
612      write(lunout,*)
613!
614      RETURN
615      END
616!
617! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
618!
619!
620      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,                 &
621     &                          ucov,vcov,temp,q,omega2)
622      USE dimphy
623      USE mod_grid_phy_lmdz
624      USE mod_phys_lmdz_para
625      USE iophy
626      USE phys_state_var_mod
627      USE iostart
628      USE write_field_phy
629      USE infotrac
630      use control_mod
631      USE comconst_mod, ONLY: im, jm, lllm
632      USE logic_mod, ONLY: fxyhypb, ysinus
633      USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn
634
635      IMPLICIT NONE
636!=======================================================
637! Ecriture du fichier de redemarrage sous format NetCDF
638!=======================================================
639!   Declarations:
640!   -------------
641      include "dimensions.h"
642!!#include "control.h"
643      include "netcdf.inc"
644
645!   Arguments:
646!   ----------
647      CHARACTER*(*) fichnom
648!Al1 plev tronque pour .nc mais plev(klev+1):=0
649      real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev)
650      real :: presnivs(klon,klev)
651      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
652      real :: q(klon,klev,nqtot),omega2(klon,klev)
653!      real :: ug(klev),vg(klev),fcoriolis
654      real :: phis(klon) 
655
656!   Variables locales pour NetCDF:
657!   ------------------------------
658      INTEGER iq
659      INTEGER length
660      PARAMETER (length = 100)
661      REAL tab_cntrl(length) ! tableau des parametres du run
662      character*4 nmq(nqtot)
663      character*12 modname
664      character*80 abort_message
665      LOGICAL found
666
667      modname = 'dyn1deta0 : '
668!!      nmq(1)="vap"
669!!      nmq(2)="cond"
670!!      do iq=3,nqtot
671!!        write(nmq(iq),'("tra",i1)') iq-2
672!!      enddo
673      DO iq = 1,nqtot
674        nmq(iq) = trim(tname(iq))
675      ENDDO
676      print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
677      CALL open_startphy(fichnom)
678      print*,'after open startphy ',fichnom,nmq
679
680!
681! Lecture des parametres de controle:
682!
683      CALL get_var("controle",tab_cntrl)
684       
685
686      im         = tab_cntrl(1)
687      jm         = tab_cntrl(2)
688      lllm       = tab_cntrl(3)
689      day_ref    = tab_cntrl(4)
690      annee_ref  = tab_cntrl(5)
691!      rad        = tab_cntrl(6)
692!      omeg       = tab_cntrl(7)
693!      g          = tab_cntrl(8)
694!      cpp        = tab_cntrl(9)
695!      kappa      = tab_cntrl(10)
696!      daysec     = tab_cntrl(11)
697!      dtvr       = tab_cntrl(12)
698!      etot0      = tab_cntrl(13)
699!      ptot0      = tab_cntrl(14)
700!      ztot0      = tab_cntrl(15)
701!      stot0      = tab_cntrl(16)
702!      ang0       = tab_cntrl(17)
703!      pa         = tab_cntrl(18)
704!      preff      = tab_cntrl(19)
705!
706!      clon       = tab_cntrl(20)
707!      clat       = tab_cntrl(21)
708!      grossismx  = tab_cntrl(22)
709!      grossismy  = tab_cntrl(23)
710!
711      IF ( tab_cntrl(24).EQ.1. )  THEN
712        fxyhypb  =.true.
713!        dzoomx   = tab_cntrl(25)
714!        dzoomy   = tab_cntrl(26)
715!        taux     = tab_cntrl(28)
716!        tauy     = tab_cntrl(29)
717      ELSE
718        fxyhypb = .false.
719        ysinus  = .false.
720        IF( tab_cntrl(27).EQ.1. ) ysinus =.true. 
721      ENDIF
722
723      day_ini = tab_cntrl(30)
724      itau_dyn = tab_cntrl(31)
725!   .................................................................
726!
727!
728!      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
729!Al1
730       Print*,'day_ref,annee_ref,day_ini,itau_dyn',                         &
731     &              day_ref,annee_ref,day_ini,itau_dyn
732
733!  Lecture des champs
734!
735      CALL get_field("play",play,found)
736      IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
737      CALL get_field("phi",phi,found)
738      IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
739      CALL get_field("phis",phis,found)
740      IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
741      CALL get_field("presnivs",presnivs,found)
742      IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
743      CALL get_field("ucov",ucov,found)
744      IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
745      CALL get_field("vcov",vcov,found)
746      IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
747      CALL get_field("temp",temp,found)
748      IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
749      CALL get_field("omega2",omega2,found)
750      IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
751      plev(1,klev+1)=0.
752      CALL get_field("plev",plev(:,1:klev),found)
753      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
754
755      Do iq=1,nqtot
756        CALL get_field("q"//nmq(iq),q(:,:,iq),found)
757        IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
758      EndDo
759
760      CALL close_startphy
761      print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
762!
763      RETURN
764      END
765!
766! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
767!
768!
769      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,           &
770     &                          ucov,vcov,temp,q,omega2)
771      USE dimphy
772      USE mod_grid_phy_lmdz
773      USE mod_phys_lmdz_para
774      USE phys_state_var_mod
775      USE iostart
776      USE infotrac
777      use control_mod
778      USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
779      USE logic_mod, ONLY: fxyhypb, ysinus
780      USE temps_mod, ONLY: annee_ref,day_end,day_ref,itau_dyn,itaufin
781
782      IMPLICIT NONE
783!=======================================================
784! Ecriture du fichier de redemarrage sous format NetCDF
785!=======================================================
786!   Declarations:
787!   -------------
788      include "dimensions.h"
789!!#include "control.h"
790      include "netcdf.inc"
791
792!   Arguments:
793!   ----------
794      CHARACTER*(*) fichnom
795!Al1 plev tronque pour .nc mais plev(klev+1):=0
796      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
797      real :: presnivs(klon,klev)
798      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
799      real :: q(klon,klev,nqtot)
800      real :: omega2(klon,klev),rho(klon,klev+1)
801!      real :: ug(klev),vg(klev),fcoriolis
802      real :: phis(klon) 
803
804!   Variables locales pour NetCDF:
805!   ------------------------------
806      INTEGER nid
807      INTEGER ierr
808      INTEGER iq,l
809      INTEGER length
810      PARAMETER (length = 100)
811      REAL tab_cntrl(length) ! tableau des parametres du run
812      character*4 nmq(nqtot)
813      character*20 modname
814      character*80 abort_message
815!
816      INTEGER nb
817      SAVE nb
818      DATA nb / 0 /
819
820      CALL open_restartphy(fichnom)
821      print*,'redm1 ',fichnom,klon,klev,nqtot
822!!      nmq(1)="vap"
823!!      nmq(2)="cond"
824!!      nmq(3)="tra1"
825!!      nmq(4)="tra2"
826      DO iq = 1,nqtot
827        nmq(iq) = trim(tname(iq))
828      ENDDO
829
830      modname = 'dyn1dredem'
831      ierr = NF_OPEN(fichnom, NF_WRITE, nid)
832      IF (ierr .NE. NF_NOERR) THEN
833         abort_message="Pb. d ouverture "//fichnom
834         CALL abort_gcm('Modele 1D',abort_message,1)
835      ENDIF
836
837      DO l=1,length
838       tab_cntrl(l) = 0.
839      ENDDO
840       tab_cntrl(1)  = FLOAT(iim)
841       tab_cntrl(2)  = FLOAT(jjm)
842       tab_cntrl(3)  = FLOAT(llm)
843       tab_cntrl(4)  = FLOAT(day_ref)
844       tab_cntrl(5)  = FLOAT(annee_ref)
845       tab_cntrl(6)  = rad
846       tab_cntrl(7)  = omeg
847       tab_cntrl(8)  = g
848       tab_cntrl(9)  = cpp
849       tab_cntrl(10) = kappa
850       tab_cntrl(11) = daysec
851       tab_cntrl(12) = dtvr
852!       tab_cntrl(13) = etot0
853!       tab_cntrl(14) = ptot0
854!       tab_cntrl(15) = ztot0
855!       tab_cntrl(16) = stot0
856!       tab_cntrl(17) = ang0
857!       tab_cntrl(18) = pa
858!       tab_cntrl(19) = preff
859!
860!    .....    parametres  pour le zoom      ......   
861
862!       tab_cntrl(20)  = clon
863!       tab_cntrl(21)  = clat
864!       tab_cntrl(22)  = grossismx
865!       tab_cntrl(23)  = grossismy
866!
867      IF ( fxyhypb )   THEN
868       tab_cntrl(24) = 1.
869!       tab_cntrl(25) = dzoomx
870!       tab_cntrl(26) = dzoomy
871       tab_cntrl(27) = 0.
872!       tab_cntrl(28) = taux
873!       tab_cntrl(29) = tauy
874      ELSE
875       tab_cntrl(24) = 0.
876!       tab_cntrl(25) = dzoomx
877!       tab_cntrl(26) = dzoomy
878       tab_cntrl(27) = 0.
879       tab_cntrl(28) = 0.
880       tab_cntrl(29) = 0.
881       IF( ysinus )  tab_cntrl(27) = 1.
882      ENDIF
883!Al1 iday_end -> day_end
884       tab_cntrl(30) = FLOAT(day_end)
885       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
886!
887      CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
888!
889
890!  Ecriture/extension de la coordonnee temps
891
892      nb = nb + 1
893
894!  Ecriture des champs
895!
896      CALL put_field("plev","p interfaces sauf la nulle",plev)
897      CALL put_field("play","",play)
898      CALL put_field("phi","geopotentielle",phi)
899      CALL put_field("phis","geopotentiell de surface",phis)
900      CALL put_field("presnivs","",presnivs)
901      CALL put_field("ucov","",ucov)
902      CALL put_field("vcov","",vcov)
903      CALL put_field("temp","",temp)
904      CALL put_field("omega2","",omega2)
905
906      Do iq=1,nqtot
907        CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",           &
908     &                                                      q(:,:,iq))
909      EndDo
910      CALL close_restartphy
911
912!
913      RETURN
914      END
915      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
916      IMPLICIT NONE
917!=======================================================================
918!   passage d'un champ de la grille scalaire a la grille physique
919!=======================================================================
920 
921!-----------------------------------------------------------------------
922!   declarations:
923!   -------------
924 
925      INTEGER im,jm,ngrid,nfield
926      REAL pdyn(im,jm,nfield)
927      REAL pfi(ngrid,nfield)
928 
929      INTEGER i,j,ifield,ig
930 
931!-----------------------------------------------------------------------
932!   calcul:
933!   -------
934 
935      DO ifield=1,nfield
936!   traitement des poles
937         DO i=1,im
938            pdyn(i,1,ifield)=pfi(1,ifield)
939            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
940         ENDDO
941 
942!   traitement des point normaux
943         DO j=2,jm-1
944            ig=2+(j-2)*(im-1)
945            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
946            pdyn(im,j,ifield)=pdyn(1,j,ifield)
947         ENDDO
948      ENDDO
949 
950      RETURN
951      END
952 
953 
954
955      SUBROUTINE abort_gcm(modname, message, ierr)
956 
957      USE IOIPSL
958!
959! Stops the simulation cleanly, closing files and printing various
960! comments
961!
962!  Input: modname = name of calling program
963!         message = stuff to print
964!         ierr    = severity of situation ( = 0 normal )
965 
966      character(len=*) modname
967      integer ierr
968      character(len=*) message
969 
970      write(*,*) 'in abort_gcm'
971      call histclo
972!     call histclo(2)
973!     call histclo(3)
974!     call histclo(4)
975!     call histclo(5)
976      write(*,*) 'out of histclo'
977      write(*,*) 'Stopping in ', modname
978      write(*,*) 'Reason = ',message
979      call getin_dump
980!
981      if (ierr .eq. 0) then
982        write(*,*) 'Everything is cool'
983      else
984        write(*,*) 'Houston, we have a problem ', ierr
985      endif
986      STOP
987      END
988      REAL FUNCTION fq_sat(kelvin, millibar)
989!
990      IMPLICIT none
991!======================================================================
992! Autheur(s): Z.X. Li (LMD/CNRS)
993! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
994!======================================================================
995! Arguments:
996! kelvin---input-R: temperature en Kelvin
997! millibar--input-R: pression en mb
998!
999! fq_sat----output-R: vapeur d'eau saturante en kg/kg
1000!======================================================================
1001!
1002      REAL kelvin, millibar
1003!
1004      REAL r2es
1005      PARAMETER (r2es=611.14 *18.0153/28.9644)
1006!
1007      REAL r3les, r3ies, r3es
1008      PARAMETER (R3LES=17.269)
1009      PARAMETER (R3IES=21.875)
1010!
1011      REAL r4les, r4ies, r4es
1012      PARAMETER (R4LES=35.86)
1013      PARAMETER (R4IES=7.66)
1014!
1015      REAL rtt
1016      PARAMETER (rtt=273.16)
1017!
1018      REAL retv
1019      PARAMETER (retv=28.9644/18.0153 - 1.0)
1020!
1021      REAL zqsat
1022      REAL temp, pres
1023!     ------------------------------------------------------------------
1024!
1025!
1026      temp = kelvin
1027      pres = millibar * 100.0
1028!      write(*,*)'kelvin,millibar=',kelvin,millibar
1029!      write(*,*)'temp,pres=',temp,pres
1030!
1031      IF (temp .LE. rtt) THEN
1032         r3es = r3ies
1033         r4es = r4ies
1034      ELSE
1035         r3es = r3les
1036         r4es = r4les
1037      ENDIF
1038!
1039      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
1040      zqsat=MIN(0.5,ZQSAT)
1041      zqsat=zqsat/(1.-retv  *zqsat)
1042!
1043      fq_sat = zqsat
1044!
1045      RETURN
1046      END
1047 
1048      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
1049      IMPLICIT NONE
1050!=======================================================================
1051!   passage d'un champ de la grille scalaire a la grille physique
1052!=======================================================================
1053 
1054!-----------------------------------------------------------------------
1055!   declarations:
1056!   -------------
1057 
1058      INTEGER im,jm,ngrid,nfield
1059      REAL pdyn(im,jm,nfield)
1060      REAL pfi(ngrid,nfield)
1061 
1062      INTEGER j,ifield,ig
1063 
1064!-----------------------------------------------------------------------
1065!   calcul:
1066!   -------
1067 
1068      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)                          &
1069     &    STOP 'probleme de dim'
1070!   traitement des poles
1071      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
1072      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
1073 
1074!   traitement des point normaux
1075      DO ifield=1,nfield
1076         DO j=2,jm-1
1077            ig=2+(j-2)*(im-1)
1078            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
1079         ENDDO
1080      ENDDO
1081 
1082      RETURN
1083      END
1084 
1085      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
1086 
1087!    Ancienne version disvert dont on a modifie nom pour utiliser
1088!    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
1089!    (MPL 18092012)
1090!
1091!    Auteur :  P. Le Van .
1092!
1093      IMPLICIT NONE
1094 
1095      include "dimensions.h"
1096      include "paramet.h"
1097!
1098!=======================================================================
1099!
1100!
1101!    s = sigma ** kappa   :  coordonnee  verticale
1102!    dsig(l)            : epaisseur de la couche l ds la coord.  s
1103!    sig(l)             : sigma a l'interface des couches l et l-1
1104!    ds(l)              : distance entre les couches l et l-1 en coord.s
1105!
1106!=======================================================================
1107!
1108      REAL pa,preff
1109      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
1110      REAL presnivs(llm)
1111!
1112!   declarations:
1113!   -------------
1114!
1115      REAL sig(llm+1),dsig(llm)
1116!
1117      INTEGER l
1118      REAL snorm
1119      REAL alpha,beta,gama,delta,deltaz,h
1120      INTEGER np,ierr
1121      REAL pi,x
1122 
1123!-----------------------------------------------------------------------
1124!
1125      pi=2.*ASIN(1.)
1126 
1127      OPEN(99,file='sigma.def',status='old',form='formatted',                   &
1128     &   iostat=ierr)
1129 
1130!-----------------------------------------------------------------------
1131!   cas 1 on lit les options dans sigma.def:
1132!   ----------------------------------------
1133 
1134      IF (ierr.eq.0) THEN
1135 
1136      print*,'WARNING!!! on lit les options dans sigma.def'
1137      READ(99,*) deltaz
1138      READ(99,*) h
1139      READ(99,*) beta
1140      READ(99,*) gama
1141      READ(99,*) delta
1142      READ(99,*) np
1143      CLOSE(99)
1144      alpha=deltaz/(llm*h)
1145!
1146 
1147       DO 1  l = 1, llm
1148       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*                    &
1149     &          ( (tanh(gama*l)/tanh(gama*llm))**np +                      &
1150     &            (1.-l/FLOAT(llm))*delta )
1151   1   CONTINUE
1152 
1153       sig(1)=1.
1154       DO 101 l=1,llm-1
1155          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
1156101    CONTINUE
1157       sig(llm+1)=0.
1158 
1159       DO 2  l = 1, llm
1160       dsig(l) = sig(l)-sig(l+1)
1161   2   CONTINUE
1162!
1163 
1164      ELSE
1165!-----------------------------------------------------------------------
1166!   cas 2 ancienne discretisation (LMD5...):
1167!   ----------------------------------------
1168 
1169      PRINT*,'WARNING!!! Ancienne discretisation verticale'
1170 
1171      h=7.
1172      snorm  = 0.
1173      DO l = 1, llm
1174         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
1175         dsig(l) = 1.0 + 7.0 * SIN(x)**2
1176         snorm = snorm + dsig(l)
1177      ENDDO
1178      snorm = 1./snorm
1179      DO l = 1, llm
1180         dsig(l) = dsig(l)*snorm
1181      ENDDO
1182      sig(llm+1) = 0.
1183      DO l = llm, 1, -1
1184         sig(l) = sig(l+1) + dsig(l)
1185      ENDDO
1186 
1187      ENDIF
1188 
1189 
1190      DO l=1,llm
1191        nivsigs(l) = FLOAT(l)
1192      ENDDO
1193 
1194      DO l=1,llmp1
1195        nivsig(l)= FLOAT(l)
1196      ENDDO
1197 
1198!
1199!    ....  Calculs  de ap(l) et de bp(l)  ....
1200!    .........................................
1201!
1202!
1203!   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
1204!
1205 
1206      bp(llmp1) =   0.
1207 
1208      DO l = 1, llm
1209!c
1210!cc    ap(l) = 0.
1211!cc    bp(l) = sig(l)
1212 
1213      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
1214      ap(l) = pa * ( sig(l) - bp(l) )
1215!
1216      ENDDO
1217      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
1218 
1219      PRINT *,' BP '
1220      PRINT *,  bp
1221      PRINT *,' AP '
1222      PRINT *,  ap
1223 
1224      DO l = 1, llm
1225       dpres(l) = bp(l) - bp(l+1)
1226       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
1227      ENDDO
1228 
1229      PRINT *,' PRESNIVS '
1230      PRINT *,presnivs
1231 
1232      RETURN
1233      END
1234
1235!======================================================================
1236       SUBROUTINE read_tsurf1d(knon,sst_out)
1237
1238! This subroutine specifies the surface temperature to be used in 1D simulations
1239
1240      USE dimphy, ONLY : klon
1241
1242      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
1243      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
1244
1245       INTEGER :: i
1246! COMMON defined in lmdz1d.F:
1247       real ts_cur
1248       common /sst_forcing/ts_cur
1249
1250       DO i = 1, knon
1251        sst_out(i) = ts_cur
1252       ENDDO
1253
1254      END SUBROUTINE read_tsurf1d
1255
1256!===============================================================
1257      subroutine advect_vert(llm,w,dt,q,plev)
1258!===============================================================
1259!   Schema amont pour l'advection verticale en 1D
1260!   w est la vitesse verticale dp/dt en Pa/s
1261!   Traitement en volumes finis
1262!   d / dt ( zm q ) = delta_z ( omega q )
1263!   d / dt ( zm ) = delta_z ( omega )
1264!   avec zm = delta_z ( p )
1265!   si * designe la valeur au pas de temps t+dt
1266!   zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1267!   zm*(l) -zm(l) = w(l+1) - w(l)
1268!   avec w=omega * dt
1269!---------------------------------------------------------------
1270      implicit none
1271! arguments
1272      integer llm
1273      real w(llm+1),q(llm),plev(llm+1),dt
1274
1275! local
1276      integer l
1277      real zwq(llm+1),zm(llm+1),zw(llm+1)
1278      real qold
1279
1280!---------------------------------------------------------------
1281
1282      do l=1,llm
1283         zw(l)=dt*w(l)
1284         zm(l)=plev(l)-plev(l+1)
1285         zwq(l)=q(l)*zw(l)
1286      enddo
1287      zwq(llm+1)=0.
1288      zw(llm+1)=0.
1289 
1290      do l=1,llm
1291         qold=q(l)
1292         q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1293         print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1294      enddo
1295
1296 
1297      return
1298      end
1299
1300!===============================================================
1301
1302
1303       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,              &
1304     &                q,temp,u,v,play)
1305!itlmd
1306!----------------------------------------------------------------------
1307!   Calcul de l'advection verticale (ascendance et subsidence) de
1308!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1309!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1310!   sans WTG rajouter une advection horizontale
1311!---------------------------------------------------------------------- 
1312        implicit none
1313#include "YOMCST.h"
1314!        argument
1315        integer llm
1316        real  omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1317        real  d_u_va(llm), d_v_va(llm)
1318        real  q(llm,3),temp(llm)
1319        real  u(llm),v(llm)
1320        real  play(llm)
1321! interne
1322        integer l
1323        real alpha,omgdown,omgup
1324
1325      do l= 1,llm
1326       if(l.eq.1) then
1327!si omgup pour la couche 1, alors tendance nulle
1328        omgdown=max(omega(2),0.0)
1329        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1330        d_t_va(l)= alpha*(omgdown)-omgdown*(temp(l)-temp(l+1))             &
1331     &       /(play(l)-play(l+1))
1332
1333        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))             
1334
1335        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))             
1336        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))             
1337
1338       
1339       elseif(l.eq.llm) then
1340        omgup=min(omega(l),0.0)
1341        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1342        d_t_va(l)= alpha*(omgup)-                                          &
1343
1344!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1345     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1346        d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1347        d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1348        d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1349       
1350       else
1351        omgup=min(omega(l),0.0)
1352        omgdown=max(omega(l+1),0.0)
1353        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1354        d_t_va(l)= alpha*(omgup+omgdown)-omgdown*(temp(l)-temp(l+1))       &
1355     &              /(play(l)-play(l+1))-                                  &
1356!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1357     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1358!      print*, '  ??? '
1359
1360        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))                            &
1361     &              /(play(l)-play(l+1))-                                  &
1362     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) 
1363        d_u_va(l)= -omgdown*(u(l)-u(l+1))                                  &
1364     &              /(play(l)-play(l+1))-                                  &
1365     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) 
1366        d_v_va(l)= -omgdown*(v(l)-v(l+1))                                  &
1367     &              /(play(l)-play(l+1))-                                  &
1368     &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1369       
1370      endif
1371         
1372      enddo
1373!fin itlmd
1374        return
1375        end
1376!       SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1377       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,                        &
1378     &                q,temp,u,v,play)
1379!itlmd
1380!----------------------------------------------------------------------
1381!   Calcul de l'advection verticale (ascendance et subsidence) de
1382!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1383!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1384!   sans WTG rajouter une advection horizontale
1385!---------------------------------------------------------------------- 
1386        implicit none
1387#include "YOMCST.h"
1388!        argument
1389        integer llm,nqtot
1390        real  omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1391!        real  d_u_va(llm), d_v_va(llm)
1392        real  q(llm,nqtot),temp(llm)
1393        real  u(llm),v(llm)
1394        real  play(llm)
1395        real cor(llm)
1396!        real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1397        real dph(llm),dqdp(llm),dtdp(llm)
1398! interne
1399        integer k
1400        real omdn,omup
1401
1402!        dudp=0.
1403!        dvdp=0.
1404        dqdp=0.
1405        dtdp=0.
1406!        d_u_va=0.
1407!        d_v_va=0.
1408
1409      cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
1410
1411
1412      do k=2,llm-1
1413
1414       dph  (k-1) = (play()- play(k-1  ))
1415!       dudp (k-1) = (u   ()- u   (k-1  ))/dph(k-1)
1416!       dvdp (k-1) = (v   ()- v   (k-1  ))/dph(k-1)
1417       dqdp (k-1) = (q   (k,1)- q   (k-1,1))/dph(k-1)
1418       dtdp (k-1) = (temp()- temp(k-1  ))/dph(k-1)
1419
1420      enddo
1421
1422!      dudp (  llm  ) = dudp ( llm-1 )
1423!      dvdp (  llm  ) = dvdp ( llm-1 )
1424      dqdp (  llm  ) = dqdp ( llm-1 )
1425      dtdp (  llm  ) = dtdp ( llm-1 )
1426
1427      do k=2,llm-1
1428      omdn=max(0.0,omega(k+1))
1429      omup=min(0.0,omega( k ))
1430
1431!      d_u_va(k)  = -omdn*dudp(k)-omup*dudp(k-1)
1432!      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
1433      d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1)
1434      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)+(omup+omdn)*cor(k)
1435      enddo
1436
1437      omdn=max(0.0,omega( 2 ))
1438      omup=min(0.0,omega(llm))
1439!      d_u_va( 1 )   = -omdn*dudp( 1 )
1440!      d_u_va(llm)   = -omup*dudp(llm)
1441!      d_v_va( 1 )   = -omdn*dvdp( 1 )
1442!      d_v_va(llm)   = -omup*dvdp(llm)
1443      d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1444      d_q_va(llm,1) = -omup*dqdp(llm)
1445      d_t_va( 1 )   = -omdn*dtdp( 1 )+omdn*cor( 1 )
1446      d_t_va(llm)   = -omup*dtdp(llm)!+omup*cor(llm)
1447
1448!      if(abs(rlat(1))>10.) then
1449!     Calculate the tendency due agestrophic motions
1450!      du_age = fcoriolis*(v-vg)
1451!      dv_age = fcoriolis*(ug-u)
1452!      endif
1453
1454!       call writefield_phy('d_t_va',d_t_va,llm)
1455
1456          return
1457         end
1458
1459!======================================================================
1460      SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga                     &
1461     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga        &
1462     &             ,ht_toga,vt_toga,hq_toga,vq_toga)
1463      implicit none
1464
1465!-------------------------------------------------------------------------
1466! Read TOGA-COARE forcing data
1467!-------------------------------------------------------------------------
1468
1469      integer nlev_toga,nt_toga
1470      real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1471      real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1472      real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1473      real w_toga(nlev_toga,nt_toga)
1474      real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1475      real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1476      character*80 fich_toga
1477
1478      integer k,ip
1479      real bid
1480
1481      integer iy,im,id,ih
1482     
1483       real plev_min
1484
1485       plev_min = 55.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1486
1487      open(21,file=trim(fich_toga),form='formatted')
1488      read(21,'(a)') 
1489      do ip = 1, nt_toga
1490      read(21,'(a)') 
1491      read(21,'(a)') 
1492      read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1493      read(21,'(a)') 
1494      read(21,'(a)') 
1495
1496       do k = 1, nlev_toga
1497         read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)          &
1498     &       ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)                     &
1499     &       ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1500
1501! conversion in SI units:
1502         t_toga(k,ip)=t_toga(k,ip)+273.15     ! K
1503         q_toga(k,ip)=q_toga(k,ip)*0.001      ! kg/kg
1504         w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1505! no water vapour tendency above 55 hPa
1506         if (plev_toga(k,ip) .lt. plev_min) then
1507          q_toga(k,ip) = 0.
1508          hq_toga(k,ip) = 0.
1509          vq_toga(k,ip) =0.
1510         endif
1511       enddo
1512
1513         ts_toga(ip)=ts_toga(ip)+273.15       ! K
1514       enddo
1515       close(21)
1516
1517  223 format(4i3,6f8.2)
1518  230 format(6f9.3,4e11.3)
1519
1520          return
1521          end
1522
1523!-------------------------------------------------------------------------
1524      SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
1525      implicit none
1526
1527!-------------------------------------------------------------------------
1528! Read I.SANDU case forcing data
1529!-------------------------------------------------------------------------
1530
1531      integer nlev_sandu,nt_sandu
1532      real ts_sandu(nt_sandu)
1533      character*80 fich_sandu
1534
1535      integer ip
1536      integer iy,im,id,ih
1537
1538      real plev_min
1539
1540      print*,'nlev_sandu',nlev_sandu
1541      plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1542
1543      open(21,file=trim(fich_sandu),form='formatted')
1544      read(21,'(a)')
1545      do ip = 1, nt_sandu
1546      read(21,'(a)')
1547      read(21,'(a)')
1548      read(21,223) iy, im, id, ih, ts_sandu(ip)
1549      print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
1550      enddo
1551      close(21)
1552
1553  223 format(4i3,f8.2)
1554
1555          return
1556          end
1557
1558!=====================================================================
1559!-------------------------------------------------------------------------
1560      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,      &
1561     & ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
1562      implicit none
1563
1564!-------------------------------------------------------------------------
1565! Read Astex case forcing data
1566!-------------------------------------------------------------------------
1567
1568      integer nlev_astex,nt_astex
1569      real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
1570      real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
1571      character*80 fich_astex
1572
1573      integer ip
1574      integer iy,im,id,ih
1575
1576       real plev_min
1577
1578      print*,'nlev_astex',nlev_astex
1579       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1580
1581      open(21,file=trim(fich_astex),form='formatted')
1582      read(21,'(a)')
1583      read(21,'(a)')
1584      do ip = 1, nt_astex
1585      read(21,'(a)')
1586      read(21,'(a)')
1587      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),             &
1588     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
1589      ts_astex(ip)=ts_astex(ip)+273.15
1590      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),             &
1591     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
1592      enddo
1593      close(21)
1594
1595  223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
1596
1597          return
1598          end
1599!=====================================================================
1600      subroutine read_twpice(fich_twpice,nlevel,ntime                       &
1601     &     ,T_srf,plev,T,q,u,v,omega                                       &
1602     &     ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1603
1604!program reading forcings of the TWP-ICE experiment
1605
1606!      use netcdf
1607
1608      implicit none
1609
1610#include "netcdf.inc"
1611
1612      integer ntime,nlevel
1613      integer l,k
1614      character*80 :: fich_twpice
1615      real*8 time(ntime)
1616      real*8 lat, lon, alt, phis
1617      real*8 lev(nlevel)
1618      real*8 plev(nlevel,ntime)
1619
1620      real*8 T(nlevel,ntime)
1621      real*8 q(nlevel,ntime),u(nlevel,ntime)
1622      real*8 v(nlevel,ntime)
1623      real*8 omega(nlevel,ntime), div(nlevel,ntime)
1624      real*8 T_adv_h(nlevel,ntime)
1625      real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1626      real*8 q_adv_v(nlevel,ntime)
1627      real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1628      real*8 s_adv_v(nlevel,ntime)
1629      real*8 p_srf_aver(ntime), p_srf_center(ntime)
1630      real*8 T_srf(ntime)
1631
1632      integer nid, ierr
1633      integer nbvar3d
1634      parameter(nbvar3d=20)
1635      integer var3didin(nbvar3d)
1636
1637      ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1638      if (ierr.NE.NF_NOERR) then
1639         write(*,*) 'ERROR: Pb opening forcings cdf file '
1640         write(*,*) NF_STRERROR(ierr)
1641         stop ""
1642      endif
1643
1644      ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1645         if(ierr/=NF_NOERR) then
1646           write(*,*) NF_STRERROR(ierr)
1647           stop 'lat'
1648         endif
1649     
1650       ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1651         if(ierr/=NF_NOERR) then
1652           write(*,*) NF_STRERROR(ierr)
1653           stop 'lon'
1654         endif
1655
1656       ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1657         if(ierr/=NF_NOERR) then
1658           write(*,*) NF_STRERROR(ierr)
1659           stop 'alt'
1660         endif
1661
1662      ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1663         if(ierr/=NF_NOERR) then
1664           write(*,*) NF_STRERROR(ierr)
1665           stop 'phis'
1666         endif
1667
1668      ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1669         if(ierr/=NF_NOERR) then
1670           write(*,*) NF_STRERROR(ierr)
1671           stop 'T'
1672         endif
1673
1674      ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1675         if(ierr/=NF_NOERR) then
1676           write(*,*) NF_STRERROR(ierr)
1677           stop 'q'
1678         endif
1679
1680      ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1681         if(ierr/=NF_NOERR) then
1682           write(*,*) NF_STRERROR(ierr)
1683           stop 'u'
1684         endif
1685
1686      ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1687         if(ierr/=NF_NOERR) then
1688           write(*,*) NF_STRERROR(ierr)
1689           stop 'v'
1690         endif
1691
1692      ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1693         if(ierr/=NF_NOERR) then
1694           write(*,*) NF_STRERROR(ierr)
1695           stop 'omega'
1696         endif
1697
1698      ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1699         if(ierr/=NF_NOERR) then
1700           write(*,*) NF_STRERROR(ierr)
1701           stop 'div'
1702         endif
1703
1704      ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1705         if(ierr/=NF_NOERR) then
1706           write(*,*) NF_STRERROR(ierr)
1707           stop 'T_adv_h'
1708         endif
1709
1710      ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1711         if(ierr/=NF_NOERR) then
1712           write(*,*) NF_STRERROR(ierr)
1713           stop 'T_adv_v'
1714         endif
1715
1716      ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1717         if(ierr/=NF_NOERR) then
1718           write(*,*) NF_STRERROR(ierr)
1719           stop 'q_adv_h'
1720         endif
1721
1722      ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1723         if(ierr/=NF_NOERR) then
1724           write(*,*) NF_STRERROR(ierr)
1725           stop 'q_adv_v'
1726         endif
1727
1728      ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1729         if(ierr/=NF_NOERR) then
1730           write(*,*) NF_STRERROR(ierr)
1731           stop 's'
1732         endif
1733
1734      ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1735         if(ierr/=NF_NOERR) then
1736           write(*,*) NF_STRERROR(ierr)
1737           stop 's_adv_h'
1738         endif
1739   
1740      ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1741         if(ierr/=NF_NOERR) then
1742           write(*,*) NF_STRERROR(ierr)
1743           stop 's_adv_v'
1744         endif
1745
1746      ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1747         if(ierr/=NF_NOERR) then
1748           write(*,*) NF_STRERROR(ierr)
1749           stop 'p_srf_aver'
1750         endif
1751
1752      ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1753         if(ierr/=NF_NOERR) then
1754           write(*,*) NF_STRERROR(ierr)
1755           stop 'p_srf_center'
1756         endif
1757
1758      ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1759         if(ierr/=NF_NOERR) then
1760           write(*,*) NF_STRERROR(ierr)
1761           stop 'T_srf'
1762         endif
1763
1764!dimensions lecture
1765      call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1766
1767!pressure
1768       do l=1,ntime
1769       do k=1,nlevel
1770          plev(k,l)=lev(k)
1771       enddo
1772       enddo
1773         
1774#ifdef NC_DOUBLE
1775         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1776#else
1777         ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1778#endif
1779         if(ierr/=NF_NOERR) then
1780            write(*,*) NF_STRERROR(ierr)
1781            stop "getvarup"
1782         endif
1783!         write(*,*)'lecture lat ok',lat
1784
1785#ifdef NC_DOUBLE
1786         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1787#else
1788         ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1789#endif
1790         if(ierr/=NF_NOERR) then
1791            write(*,*) NF_STRERROR(ierr)
1792            stop "getvarup"
1793         endif
1794!         write(*,*)'lecture lon ok',lon
1795 
1796#ifdef NC_DOUBLE
1797         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1798#else
1799         ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1800#endif
1801         if(ierr/=NF_NOERR) then
1802            write(*,*) NF_STRERROR(ierr)
1803            stop "getvarup"
1804         endif
1805!          write(*,*)'lecture alt ok',alt
1806 
1807#ifdef NC_DOUBLE
1808         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1809#else
1810         ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1811#endif
1812         if(ierr/=NF_NOERR) then
1813            write(*,*) NF_STRERROR(ierr)
1814            stop "getvarup"
1815         endif
1816!          write(*,*)'lecture phis ok',phis
1817         
1818#ifdef NC_DOUBLE
1819         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1820#else
1821         ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1822#endif
1823         if(ierr/=NF_NOERR) then
1824            write(*,*) NF_STRERROR(ierr)
1825            stop "getvarup"
1826         endif
1827!         write(*,*)'lecture T ok'
1828
1829#ifdef NC_DOUBLE
1830         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1831#else
1832         ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1833#endif
1834         if(ierr/=NF_NOERR) then
1835            write(*,*) NF_STRERROR(ierr)
1836            stop "getvarup"
1837         endif
1838!         write(*,*)'lecture q ok'
1839!q in kg/kg
1840       do l=1,ntime
1841       do k=1,nlevel
1842          q(k,l)=q(k,l)/1000.
1843       enddo
1844       enddo
1845#ifdef NC_DOUBLE
1846         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1847#else
1848         ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1849#endif
1850         if(ierr/=NF_NOERR) then
1851            write(*,*) NF_STRERROR(ierr)
1852            stop "getvarup"
1853         endif
1854!         write(*,*)'lecture u ok'
1855
1856#ifdef NC_DOUBLE
1857         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1858#else
1859         ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1860#endif
1861         if(ierr/=NF_NOERR) then
1862            write(*,*) NF_STRERROR(ierr)
1863            stop "getvarup"
1864         endif
1865!         write(*,*)'lecture v ok'
1866
1867#ifdef NC_DOUBLE
1868         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1869#else
1870         ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1871#endif
1872         if(ierr/=NF_NOERR) then
1873            write(*,*) NF_STRERROR(ierr)
1874            stop "getvarup"
1875         endif
1876!         write(*,*)'lecture omega ok'
1877!omega in mb/hour
1878       do l=1,ntime
1879       do k=1,nlevel
1880          omega(k,l)=omega(k,l)*100./3600.
1881       enddo
1882       enddo
1883
1884#ifdef NC_DOUBLE
1885         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1886#else
1887         ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1888#endif
1889         if(ierr/=NF_NOERR) then
1890            write(*,*) NF_STRERROR(ierr)
1891            stop "getvarup"
1892         endif
1893!         write(*,*)'lecture div ok'
1894
1895#ifdef NC_DOUBLE
1896         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1897#else
1898         ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1899#endif
1900         if(ierr/=NF_NOERR) then
1901            write(*,*) NF_STRERROR(ierr)
1902            stop "getvarup"
1903         endif
1904!         write(*,*)'lecture T_adv_h ok'
1905!T adv in K/s
1906       do l=1,ntime
1907       do k=1,nlevel
1908          T_adv_h(k,l)=T_adv_h(k,l)/3600.
1909       enddo
1910       enddo
1911
1912
1913#ifdef NC_DOUBLE
1914         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1915#else
1916         ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1917#endif
1918         if(ierr/=NF_NOERR) then
1919            write(*,*) NF_STRERROR(ierr)
1920            stop "getvarup"
1921         endif
1922!         write(*,*)'lecture T_adv_v ok'
1923!T adv in K/s
1924       do l=1,ntime
1925       do k=1,nlevel
1926          T_adv_v(k,l)=T_adv_v(k,l)/3600.
1927       enddo
1928       enddo
1929
1930#ifdef NC_DOUBLE
1931         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1932#else
1933         ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1934#endif
1935         if(ierr/=NF_NOERR) then
1936            write(*,*) NF_STRERROR(ierr)
1937            stop "getvarup"
1938         endif
1939!         write(*,*)'lecture q_adv_h ok'
1940!q adv in kg/kg/s
1941       do l=1,ntime
1942       do k=1,nlevel
1943          q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1944       enddo
1945       enddo
1946
1947
1948#ifdef NC_DOUBLE
1949         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1950#else
1951         ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1952#endif
1953         if(ierr/=NF_NOERR) then
1954            write(*,*) NF_STRERROR(ierr)
1955            stop "getvarup"
1956         endif
1957!         write(*,*)'lecture q_adv_v ok'
1958!q adv in kg/kg/s
1959       do l=1,ntime
1960       do k=1,nlevel
1961          q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1962       enddo
1963       enddo
1964
1965
1966#ifdef NC_DOUBLE
1967         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1968#else
1969         ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1970#endif
1971         if(ierr/=NF_NOERR) then
1972            write(*,*) NF_STRERROR(ierr)
1973            stop "getvarup"
1974         endif
1975
1976#ifdef NC_DOUBLE
1977         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1978#else
1979         ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1980#endif
1981         if(ierr/=NF_NOERR) then
1982            write(*,*) NF_STRERROR(ierr)
1983            stop "getvarup"
1984         endif
1985
1986#ifdef NC_DOUBLE
1987         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1988#else
1989         ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1990#endif
1991         if(ierr/=NF_NOERR) then
1992            write(*,*) NF_STRERROR(ierr)
1993            stop "getvarup"
1994         endif
1995
1996#ifdef NC_DOUBLE
1997         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1998#else
1999         ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
2000#endif
2001         if(ierr/=NF_NOERR) then
2002            write(*,*) NF_STRERROR(ierr)
2003            stop "getvarup"
2004         endif
2005
2006#ifdef NC_DOUBLE
2007         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
2008#else
2009         ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
2010#endif
2011         if(ierr/=NF_NOERR) then
2012            write(*,*) NF_STRERROR(ierr)
2013            stop "getvarup"
2014         endif
2015
2016#ifdef NC_DOUBLE
2017         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
2018#else
2019         ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
2020#endif
2021         if(ierr/=NF_NOERR) then
2022            write(*,*) NF_STRERROR(ierr)
2023            stop "getvarup"
2024         endif
2025!         write(*,*)'lecture T_srf ok', T_srf
2026
2027         return 
2028         end subroutine read_twpice
2029!=====================================================================
2030         subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
2031
2032!         use netcdf
2033
2034         implicit none
2035#include "netcdf.inc"
2036         integer nid,ttm,llm
2037         real*8 time(ttm)
2038         real*8 lev(llm)
2039         integer ierr
2040
2041         integer timevar,levvar
2042         integer timelen,levlen
2043         integer timedimin,levdimin
2044
2045! Control & lecture on dimensions
2046! ===============================
2047         ierr=NF_INQ_DIMID(nid,"time",timedimin)
2048         ierr=NF_INQ_VARID(nid,"time",timevar)
2049         if (ierr.NE.NF_NOERR) then
2050            write(*,*) 'ERROR: Field <time> is missing'
2051            stop "" 
2052         endif
2053         ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
2054
2055         ierr=NF_INQ_DIMID(nid,"lev",levdimin)
2056         ierr=NF_INQ_VARID(nid,"lev",levvar)
2057         if (ierr.NE.NF_NOERR) then
2058             write(*,*) 'ERROR: Field <lev> is lacking'
2059             stop "" 
2060         endif
2061         ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
2062
2063         if((timelen/=ttm).or.(levlen/=llm)) then
2064            write(*,*) 'ERROR: Not the good lenght for axis'
2065            write(*,*) 'longitude: ',timelen,ttm+1
2066            write(*,*) 'latitude: ',levlen,llm
2067            stop "" 
2068         endif
2069
2070!#ifdef NC_DOUBLE
2071         ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
2072         ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
2073!#else
2074!        ierr = NF_GET_VAR_REAL(nid,timevar,time)
2075!        ierr = NF_GET_VAR_REAL(nid,levvar,lev)
2076!#endif
2077
2078       return
2079       end
2080!=====================================================================
2081
2082       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof          &
2083     &         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof                &
2084     &         ,omega_prof,o3mmr_prof                                      &
2085     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                      &
2086     &         ,omega_mod,o3mmr_mod,mxcalc)
2087
2088       implicit none
2089
2090#include "dimensions.h"
2091
2092!-------------------------------------------------------------------------
2093! Vertical interpolation of SANDUREF forcing data onto model levels
2094!-------------------------------------------------------------------------
2095
2096       integer nlevmax
2097       parameter (nlevmax=41)
2098       integer nlev_sandu,mxcalc
2099!       real play(llm), plev_prof(nlevmax)
2100!       real t_prof(nlevmax),q_prof(nlevmax)
2101!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2102!       real ht_prof(nlevmax),vt_prof(nlevmax)
2103!       real hq_prof(nlevmax),vq_prof(nlevmax)
2104
2105       real play(llm), plev_prof(nlev_sandu)
2106       real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
2107       real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
2108       real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
2109
2110       real t_mod(llm),thl_mod(llm),q_mod(llm)
2111       real u_mod(llm),v_mod(llm), w_mod(llm)
2112       real omega_mod(llm),o3mmr_mod(llm)
2113
2114       integer l,k,k1,k2
2115       real frac,frac1,frac2,fact
2116
2117       do l = 1, llm
2118
2119        if (play(l).ge.plev_prof(nlev_sandu)) then
2120
2121        mxcalc=l
2122         k1=0
2123         k2=0
2124
2125         if (play(l).le.plev_prof(1)) then
2126
2127         do k = 1, nlev_sandu-1
2128          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2129            k1=k
2130            k2=k+1
2131          endif
2132         enddo
2133
2134         if (k1.eq.0 .or. k2.eq.0) then
2135          write(*,*) 'PB! k1, k2 = ',k1,k2
2136          write(*,*) 'l,play(l) = ',l,play(l)/100
2137         do k = 1, nlev_sandu-1
2138          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2139         enddo
2140         endif
2141
2142         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2143         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2144         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
2145         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2146         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2147         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2148         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2149         omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
2150         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
2151
2152         else !play>plev_prof(1)
2153
2154         k1=1
2155         k2=2
2156         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2157         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2158         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
2159         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
2160         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
2161         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2162         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2163         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2164         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
2165         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
2166
2167         endif ! play.le.plev_prof(1)
2168
2169        else ! above max altitude of forcing file
2170
2171!jyg
2172         fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
2173         fact = max(fact,0.)                                           !jyg
2174         fact = exp(-fact)                                             !jyg
2175         t_mod(l)= t_prof(nlev_sandu)                                   !jyg
2176         thl_mod(l)= thl_prof(nlev_sandu)                                   !jyg
2177         q_mod(l)= q_prof(nlev_sandu)*fact                              !jyg
2178         u_mod(l)= u_prof(nlev_sandu)*fact                              !jyg
2179         v_mod(l)= v_prof(nlev_sandu)*fact                              !jyg
2180         w_mod(l)= w_prof(nlev_sandu)*fact                              !jyg
2181         omega_mod(l)= omega_prof(nlev_sandu)*fact                      !jyg
2182         o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact                      !jyg
2183
2184        endif ! play
2185
2186       enddo ! l
2187
2188       do l = 1,llm
2189!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
2190!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
2191       enddo
2192
2193          return
2194          end
2195!=====================================================================
2196       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof          &
2197     &         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof      &
2198     &         ,w_prof,tke_prof,o3mmr_prof                                 &
2199     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod       &
2200     &         ,tke_mod,o3mmr_mod,mxcalc)
2201
2202       implicit none
2203
2204#include "dimensions.h"
2205
2206!-------------------------------------------------------------------------
2207! Vertical interpolation of Astex forcing data onto model levels
2208!-------------------------------------------------------------------------
2209
2210       integer nlevmax
2211       parameter (nlevmax=41)
2212       integer nlev_astex,mxcalc
2213!       real play(llm), plev_prof(nlevmax)
2214!       real t_prof(nlevmax),qv_prof(nlevmax)
2215!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2216!       real ht_prof(nlevmax),vt_prof(nlevmax)
2217!       real hq_prof(nlevmax),vq_prof(nlevmax)
2218
2219       real play(llm), plev_prof(nlev_astex)
2220       real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
2221       real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
2222       real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
2223       real qt_prof(nlev_astex),tke_prof(nlev_astex)
2224
2225       real t_mod(llm),thl_mod(llm),qv_mod(llm)
2226       real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
2227       real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
2228
2229       integer l,k,k1,k2
2230       real frac,frac1,frac2,fact
2231
2232       do l = 1, llm
2233
2234        if (play(l).ge.plev_prof(nlev_astex)) then
2235
2236        mxcalc=l
2237         k1=0
2238         k2=0
2239
2240         if (play(l).le.plev_prof(1)) then
2241
2242         do k = 1, nlev_astex-1
2243          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2244            k1=k
2245            k2=k+1
2246          endif
2247         enddo
2248
2249         if (k1.eq.0 .or. k2.eq.0) then
2250          write(*,*) 'PB! k1, k2 = ',k1,k2
2251          write(*,*) 'l,play(l) = ',l,play(l)/100
2252         do k = 1, nlev_astex-1
2253          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2254         enddo
2255         endif
2256
2257         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2258         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2259         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
2260         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
2261         ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
2262         qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
2263         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2264         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2265         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2266         tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
2267         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
2268
2269         else !play>plev_prof(1)
2270
2271         k1=1
2272         k2=2
2273         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2274         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2275         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
2276         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
2277         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
2278         ql_mod(l)= frac1*ql_prof(k1) - frac2*ql_prof(k2)
2279         qt_mod(l)= frac1*qt_prof(k1) - frac2*qt_prof(k2)
2280         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2281         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2282         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2283         tke_mod(l)= frac1*tke_prof(k1) - frac2*tke_prof(k2)
2284         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
2285
2286         endif ! play.le.plev_prof(1)
2287
2288        else ! above max altitude of forcing file
2289
2290!jyg
2291         fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
2292         fact = max(fact,0.)                                           !jyg
2293         fact = exp(-fact)                                             !jyg
2294         t_mod(l)= t_prof(nlev_astex)                                   !jyg
2295         thl_mod(l)= thl_prof(nlev_astex)                                   !jyg
2296         qv_mod(l)= qv_prof(nlev_astex)*fact                              !jyg
2297         ql_mod(l)= ql_prof(nlev_astex)*fact                              !jyg
2298         qt_mod(l)= qt_prof(nlev_astex)*fact                              !jyg
2299         u_mod(l)= u_prof(nlev_astex)*fact                              !jyg
2300         v_mod(l)= v_prof(nlev_astex)*fact                              !jyg
2301         w_mod(l)= w_prof(nlev_astex)*fact                              !jyg
2302         tke_mod(l)= tke_prof(nlev_astex)*fact                              !jyg
2303         o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact                      !jyg
2304
2305        endif ! play
2306
2307       enddo ! l
2308
2309       do l = 1,llm
2310!      print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
2311!    $        l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
2312       enddo
2313
2314          return
2315          end
2316
2317!======================================================================
2318      SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play                &
2319     &             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico             &
2320     &             ,dth_dyn,dqh_dyn)
2321      implicit none
2322
2323!-------------------------------------------------------------------------
2324! Read RICO forcing data
2325!-------------------------------------------------------------------------
2326#include "dimensions.h"
2327
2328
2329      integer nlev_rico
2330      real ts_rico,ps_rico
2331      real t_rico(llm),q_rico(llm)
2332      real u_rico(llm),v_rico(llm)
2333      real w_rico(llm)
2334      real dth_dyn(llm)
2335      real dqh_dyn(llm)
2336     
2337
2338      real play(llm),zlay(llm)
2339     
2340
2341      real prico(nlev_rico),zrico(nlev_rico)
2342
2343      character*80 fich_rico
2344
2345      integer k,l
2346
2347     
2348      print*,fich_rico
2349      open(21,file=trim(fich_rico),form='formatted')
2350        do k=1,llm
2351      zlay(k)=0.
2352         enddo
2353     
2354        read(21,*) ps_rico,ts_rico
2355        prico(1)=ps_rico
2356        zrico(1)=0.0
2357      do l=2,nlev_rico
2358        read(21,*) k,prico(l),zrico(l)
2359      enddo
2360       close(21)
2361
2362      do k=1,llm
2363        do l=1,80
2364          if(prico(l)>play(k)) then
2365              if(play(k)>prico(l+1)) then
2366                zlay(k)=zrico(l)+(play(k)-prico(l)) *                      &
2367     &              (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2368              else 
2369                zlay(k)=zrico(l)+(play(k)-prico(80))*                      &
2370     &              (zrico(81)-zrico(80))/(prico(81)-prico(80))
2371              endif
2372          endif
2373        enddo
2374        print*,k,zlay(k)
2375        ! U
2376        if(0 < zlay(k) .and. zlay(k) < 4000) then
2377          u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2378        elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2379       u_rico(k)=  -1.9 + (30.0 + 1.9) /                                   &
2380     &          (12000 - 4000) * (zlay(k) - 4000)
2381        elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2382          u_rico(k)=30.0
2383        elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2384          u_rico(k)=30.0 - (30.0) /                                        &
2385     & (20000 - 13000) * (zlay(k) - 13000)
2386        else
2387          u_rico(k)=0.0
2388        endif
2389
2390!Q_v
2391        if(0 < zlay(k) .and. zlay(k) < 740) then
2392          q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2393        elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2394          q_rico(k)=13.8 + (2.4 - 13.8) /                                   &
2395     &          (3260 - 740) * (zlay(k) - 740)
2396        elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2397          q_rico(k)=2.4 + (1.8 - 2.4) /                                    &
2398     &               (4000 - 3260) * (zlay(k) - 3260)
2399        elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2400          q_rico(k)=1.8 + (0 - 1.8) /                                      &
2401     &             (9000 - 4000) * (zlay(k) - 4000)
2402        else
2403          q_rico(k)=0.0
2404        endif
2405
2406!T
2407        if(0 < zlay(k) .and. zlay(k) < 740) then
2408          t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2409        elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2410          t_rico(k)=292.0 + (278.0 - 292.0) /                              &                       
2411     &       (4000 - 740) * (zlay(k) - 740)
2412        elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2413          t_rico(k)=278.0 + (203.0 - 278.0) /                              &
2414     &       (15000 - 4000) * (zlay(k) - 4000)
2415        elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2416          t_rico(k)=203.0 + (194.0 - 203.0) /                              & 
2417     &       (17500 - 15000)* (zlay(k) - 15000)
2418        elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2419          t_rico(k)=194.0 + (206.0 - 194.0) /                              &
2420     &       (20000 - 17500)* (zlay(k) - 17500)
2421        elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2422          t_rico(k)=206.0 + (270.0 - 206.0) /                              & 
2423     &        (60000 - 20000)* (zlay(k) - 20000)
2424        endif
2425
2426! W
2427        if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2428          w_rico(k)=- (0.005/2260) * zlay(k)
2429        elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2430          w_rico(k)=- 0.005
2431        elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2432       w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2433        else
2434          w_rico(k)=0.0
2435        endif
2436
2437! dThrz+dTsw0+dTlw0
2438        if(0 < zlay(k) .and. zlay(k) < 4000) then
2439          dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/                     &
2440     &               (86400*4000) * zlay(k)
2441        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2442          dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /                           &
2443     &           (86400*(5000 - 4000)) * (zlay(k) - 4000)
2444        else
2445          dth_dyn(k)=0.0
2446        endif
2447! dQhrz
2448        if(0 < zlay(k) .and. zlay(k) < 3000) then
2449          dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/                         &
2450     &                    (86400*3000) * (zlay(k))
2451        elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2452          dqh_dyn(k)=0.345 / 86400
2453        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2454          dqh_dyn(k)=0.345 / 86400 +                                       &
2455     &   (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2456        else
2457          dqh_dyn(k)=0.0
2458        endif
2459
2460!?        if(play(k)>6e4) then
2461!?          ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2462!?        elseif((play(k)>3e4).and.(play(k)<6e4)) then
2463!?          ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2464!?                          *(6e4-play(k))/(6e4-3e4)
2465!?        else
2466!?          ratqs0(1,k)=ratqshaut
2467!?        endif
2468
2469      enddo
2470
2471      do k=1,llm
2472      q_rico(k)=q_rico(k)/1e3 
2473      dqh_dyn(k)=dqh_dyn(k)/1e3
2474      v_rico(k)=-3.8
2475      enddo
2476
2477          return
2478          end
2479
2480!======================================================================
2481        SUBROUTINE interp_sandu_time(day,day1,annee_ref                    &
2482     &             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu         &
2483     &             ,nlev_sandu,ts_sandu,ts_prof)
2484        implicit none
2485
2486!---------------------------------------------------------------------------------------
2487! Time interpolation of a 2D field to the timestep corresponding to day
2488!
2489! day: current julian day (e.g. 717538.2)
2490! day1: first day of the simulation
2491! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
2492! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
2493!---------------------------------------------------------------------------------------
2494! inputs:
2495        integer annee_ref
2496        integer nt_sandu,nlev_sandu
2497        integer year_ini_sandu
2498        real day, day1,day_ini_sandu,dt_sandu
2499        real ts_sandu(nt_sandu)
2500! outputs:
2501        real ts_prof
2502! local:
2503        integer it_sandu1, it_sandu2
2504        real timeit,time_sandu1,time_sandu2,frac
2505! Check that initial day of the simulation consistent with SANDU period:
2506       if (annee_ref.ne.2006 ) then
2507        print*,'Pour SANDUREF, annee_ref doit etre 2006 '
2508        print*,'Changer annee_ref dans run.def'
2509        stop
2510       endif
2511!      if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
2512!       print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
2513!       print*,'Changer dayref dans run.def'
2514!       stop
2515!      endif
2516
2517! Determine timestep relative to the 1st day of TOGA-COARE:
2518!       timeit=(day-day1)*86400.
2519!       if (annee_ref.eq.1992) then
2520!        timeit=(day-day_ini_sandu)*86400.
2521!       else
2522!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2523!       endif
2524      timeit=(day-day_ini_sandu)*86400
2525
2526! Determine the closest observation times:
2527       it_sandu1=INT(timeit/dt_sandu)+1
2528       it_sandu2=it_sandu1 + 1
2529       time_sandu1=(it_sandu1-1)*dt_sandu
2530       time_sandu2=(it_sandu2-1)*dt_sandu
2531       print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
2532       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',              &
2533     &          it_sandu1,it_sandu2,time_sandu1,time_sandu2
2534
2535       if (it_sandu1 .ge. nt_sandu) then
2536        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '          &
2537     &        ,day,it_sandu1,it_sandu2,timeit/86400.
2538        stop
2539       endif
2540
2541! time interpolation:
2542       frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
2543       frac=max(frac,0.0)
2544
2545       ts_prof = ts_sandu(it_sandu2)                                       &
2546     &          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
2547
2548         print*,                                                           &
2549     &'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',       &
2550     &day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,                  &
2551     &it_sandu2,ts_prof
2552
2553        return
2554        END
2555!=====================================================================
2556!-------------------------------------------------------------------------
2557      SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,                &
2558     & sens,flat,adv_theta,rad_theta,adv_qt)
2559      implicit none
2560
2561!-------------------------------------------------------------------------
2562! Read ARM_CU case forcing data
2563!-------------------------------------------------------------------------
2564
2565      integer nlev_armcu,nt_armcu
2566      real sens(nt_armcu),flat(nt_armcu)
2567      real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2568      character*80 fich_armcu
2569
2570      integer ip
2571
2572      integer iy,im,id,ih,in
2573
2574      print*,'nlev_armcu',nlev_armcu
2575
2576      open(21,file=trim(fich_armcu),form='formatted')
2577      read(21,'(a)')
2578      do ip = 1, nt_armcu
2579      read(21,'(a)')
2580      read(21,'(a)')
2581      read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),                  &
2582     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2583      print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),               &
2584     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2585      enddo
2586      close(21)
2587
2588  223 format(5i3,5f8.3)
2589
2590          return
2591          end
2592
2593!=====================================================================
2594       SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof            &
2595     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                         &
2596     &         ,ht_prof,vt_prof,hq_prof,vq_prof                            &
2597     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                              &
2598     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2599 
2600       implicit none
2601 
2602#include "dimensions.h"
2603
2604!-------------------------------------------------------------------------
2605! Vertical interpolation of TOGA-COARE forcing data onto model levels
2606!-------------------------------------------------------------------------
2607 
2608       integer nlevmax
2609       parameter (nlevmax=41)
2610       integer nlev_toga,mxcalc
2611!       real play(llm), plev_prof(nlevmax) 
2612!       real t_prof(nlevmax),q_prof(nlevmax)
2613!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2614!       real ht_prof(nlevmax),vt_prof(nlevmax)
2615!       real hq_prof(nlevmax),vq_prof(nlevmax)
2616 
2617       real play(llm), plev_prof(nlev_toga) 
2618       real t_prof(nlev_toga),q_prof(nlev_toga)
2619       real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2620       real ht_prof(nlev_toga),vt_prof(nlev_toga)
2621       real hq_prof(nlev_toga),vq_prof(nlev_toga)
2622 
2623       real t_mod(llm),q_mod(llm)
2624       real u_mod(llm),v_mod(llm), w_mod(llm)
2625       real ht_mod(llm),vt_mod(llm)
2626       real hq_mod(llm),vq_mod(llm)
2627 
2628       integer l,k,k1,k2
2629       real frac,frac1,frac2,fact
2630 
2631       do l = 1, llm
2632
2633        if (play(l).ge.plev_prof(nlev_toga)) then
2634 
2635        mxcalc=l
2636         k1=0
2637         k2=0
2638
2639         if (play(l).le.plev_prof(1)) then
2640
2641         do k = 1, nlev_toga-1
2642          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2643            k1=k
2644            k2=k+1
2645          endif
2646         enddo
2647
2648         if (k1.eq.0 .or. k2.eq.0) then
2649          write(*,*) 'PB! k1, k2 = ',k1,k2
2650          write(*,*) 'l,play(l) = ',l,play(l)/100
2651         do k = 1, nlev_toga-1
2652          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2653         enddo
2654         endif
2655
2656         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2657         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2658         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2659         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2660         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2661         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2662         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2663         vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2664         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2665         vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2666     
2667         else !play>plev_prof(1)
2668
2669         k1=1
2670         k2=2
2671         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2672         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2673         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
2674         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
2675         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2676         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2677         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2678         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
2679         vt_mod(l)= frac1*vt_prof(k1) - frac2*vt_prof(k2)
2680         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
2681         vq_mod(l)= frac1*vq_prof(k1) - frac2*vq_prof(k2)
2682
2683         endif ! play.le.plev_prof(1)
2684
2685        else ! above max altitude of forcing file
2686 
2687!jyg
2688         fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2689         fact = max(fact,0.)                                           !jyg
2690         fact = exp(-fact)                                             !jyg
2691         t_mod(l)= t_prof(nlev_toga)                                   !jyg
2692         q_mod(l)= q_prof(nlev_toga)*fact                              !jyg
2693         u_mod(l)= u_prof(nlev_toga)*fact                              !jyg
2694         v_mod(l)= v_prof(nlev_toga)*fact                              !jyg
2695         w_mod(l)= 0.0                                                 !jyg
2696         ht_mod(l)= ht_prof(nlev_toga)                                 !jyg
2697         vt_mod(l)= vt_prof(nlev_toga)                                 !jyg
2698         hq_mod(l)= hq_prof(nlev_toga)*fact                            !jyg
2699         vq_mod(l)= vq_prof(nlev_toga)*fact                            !jyg
2700 
2701        endif ! play
2702 
2703       enddo ! l
2704
2705!       do l = 1,llm
2706!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2707!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2708!       enddo
2709 
2710          return
2711          end
2712 
2713!=====================================================================
2714       SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas            &
2715     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
2716     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
2717     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
2718     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
2719     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
2720     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
2721 
2722       implicit none
2723 
2724#include "dimensions.h"
2725
2726!-------------------------------------------------------------------------
2727! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
2728!-------------------------------------------------------------------------
2729 
2730       integer nlevmax
2731       parameter (nlevmax=41)
2732       integer nlev_cas,mxcalc
2733!       real play(llm), plev_prof(nlevmax) 
2734!       real t_prof(nlevmax),q_prof(nlevmax)
2735!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2736!       real ht_prof(nlevmax),vt_prof(nlevmax)
2737!       real hq_prof(nlevmax),vq_prof(nlevmax)
2738 
2739       real play(llm), plev_prof_cas(nlev_cas) 
2740       real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
2741       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
2742       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
2743       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
2744       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
2745       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
2746       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
2747 
2748       real t_mod_cas(llm),q_mod_cas(llm)
2749       real u_mod_cas(llm),v_mod_cas(llm)
2750       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
2751       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
2752       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
2753       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
2754       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
2755 
2756       integer l,k,k1,k2
2757       real frac,frac1,frac2,fact
2758 
2759       do l = 1, llm
2760
2761        if (play(l).ge.plev_prof_cas(nlev_cas)) then
2762 
2763        mxcalc=l
2764         k1=0
2765         k2=0
2766
2767         if (play(l).le.plev_prof_cas(1)) then
2768
2769         do k = 1, nlev_cas-1
2770          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
2771            k1=k
2772            k2=k+1
2773          endif
2774         enddo
2775
2776         if (k1.eq.0 .or. k2.eq.0) then
2777          write(*,*) 'PB! k1, k2 = ',k1,k2
2778          write(*,*) 'l,play(l) = ',l,play(l)/100
2779         do k = 1, nlev_cas-1
2780          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
2781         enddo
2782         endif
2783
2784         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
2785         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
2786         q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
2787         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
2788         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
2789         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
2790         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
2791         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
2792         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
2793         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
2794         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
2795         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
2796         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
2797         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
2798         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
2799         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
2800         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
2801         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
2802         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
2803         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
2804         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
2805     
2806         else !play>plev_prof_cas(1)
2807
2808         k1=1
2809         k2=2
2810         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
2811         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
2812         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
2813         q_mod_cas(l)= frac1*q_prof_cas(k1) - frac2*q_prof_cas(k2)
2814         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
2815         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
2816         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
2817         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
2818         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
2819         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
2820         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
2821         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
2822         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
2823         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
2824         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
2825         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
2826         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
2827         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
2828         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
2829         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
2830         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
2831         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
2832
2833         endif ! play.le.plev_prof_cas(1)
2834
2835        else ! above max altitude of forcing file
2836 
2837!jyg
2838         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
2839         fact = max(fact,0.)                                           !jyg
2840         fact = exp(-fact)                                             !jyg
2841         t_mod_cas(l)= t_prof_cas(nlev_cas)                                   !jyg
2842         q_mod_cas(l)= q_prof_cas(nlev_cas)*fact                              !jyg
2843         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                              !jyg
2844         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                              !jyg
2845         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                              !jyg
2846         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                              !jyg
2847         w_mod_cas(l)= 0.0                                                 !jyg
2848         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
2849         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                            !jyg
2850         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                            !jyg
2851         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
2852         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                            !jyg
2853         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                            !jyg
2854         dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
2855         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                                 !jyg
2856         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                                 !jyg
2857         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
2858         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
2859         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
2860         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact                      !jyg
2861 
2862        endif ! play
2863 
2864       enddo ! l
2865
2866!       do l = 1,llm
2867!       print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
2868!     $        l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
2869!       enddo
2870 
2871          return
2872          end
2873!***************************************************************************** 
2874!=====================================================================
2875       SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof   &
2876     &         ,th_prof,qv_prof,u_prof,v_prof,o3_prof                     &
2877     &         ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof         &
2878     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                          &
2879     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
2880 
2881       implicit none
2882 
2883#include "dimensions.h"
2884
2885!-------------------------------------------------------------------------
2886! Vertical interpolation of Dice forcing data onto model levels
2887!-------------------------------------------------------------------------
2888 
2889       integer nlevmax
2890       parameter (nlevmax=41)
2891       integer nlev_dice,mxcalc,nt_dice
2892 
2893       real play(llm), plev_prof(nlev_dice) 
2894       real th_prof(nlev_dice),qv_prof(nlev_dice)
2895       real u_prof(nlev_dice),v_prof(nlev_dice) 
2896       real o3_prof(nlev_dice)
2897       real ht_prof(nlev_dice),hq_prof(nlev_dice)
2898       real hu_prof(nlev_dice),hv_prof(nlev_dice)
2899       real w_prof(nlev_dice),omega_prof(nlev_dice)
2900 
2901       real th_mod(llm),qv_mod(llm)
2902       real u_mod(llm),v_mod(llm), o3_mod(llm)
2903       real ht_mod(llm),hq_mod(llm)
2904       real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
2905 
2906       integer l,k,k1,k2,kp
2907       real aa,frac,frac1,frac2,fact
2908 
2909       do l = 1, llm
2910
2911        if (play(l).ge.plev_prof(nlev_dice)) then
2912 
2913        mxcalc=l
2914         k1=0
2915         k2=0
2916
2917         if (play(l).le.plev_prof(1)) then
2918
2919         do k = 1, nlev_dice-1
2920          if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
2921            k1=k
2922            k2=k+1
2923          endif
2924         enddo
2925
2926         if (k1.eq.0 .or. k2.eq.0) then
2927          write(*,*) 'PB! k1, k2 = ',k1,k2
2928          write(*,*) 'l,play(l) = ',l,play(l)/100
2929         do k = 1, nlev_dice-1
2930          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2931         enddo
2932         endif
2933
2934         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2935         th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
2936         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
2937         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2938         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2939         o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
2940         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2941         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2942         hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
2943         hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
2944         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2945         omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
2946     
2947         else !play>plev_prof(1)
2948
2949         k1=1
2950         k2=2
2951         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2952         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2953         th_mod(l)= frac1*th_prof(k1) - frac2*th_prof(k2)
2954         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
2955         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2956         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2957         o3_mod(l)= frac1*o3_prof(k1) - frac2*o3_prof(k2)
2958         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
2959         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
2960         hu_mod(l)= frac1*hu_prof(k1) - frac2*hu_prof(k2)
2961         hv_mod(l)= frac1*hv_prof(k1) - frac2*hv_prof(k2)
2962         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2963         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
2964
2965         endif ! play.le.plev_prof(1)
2966
2967        else ! above max altitude of forcing file
2968 
2969!jyg
2970         fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
2971         fact = max(fact,0.)                                           !jyg
2972         fact = exp(-fact)                                             !jyg
2973         th_mod(l)= th_prof(nlev_dice)                                 !jyg
2974         qv_mod(l)= qv_prof(nlev_dice)*fact                            !jyg
2975         u_mod(l)= u_prof(nlev_dice)*fact                              !jyg
2976         v_mod(l)= v_prof(nlev_dice)*fact                              !jyg
2977         o3_mod(l)= o3_prof(nlev_dice)*fact                            !jyg
2978         ht_mod(l)= ht_prof(nlev_dice)                                 !jyg
2979         hq_mod(l)= hq_prof(nlev_dice)*fact                            !jyg
2980         hu_mod(l)= hu_prof(nlev_dice)                                 !jyg
2981         hv_mod(l)= hv_prof(nlev_dice)                                 !jyg
2982         w_mod(l)= 0.                                                  !jyg
2983         omega_mod(l)= 0.                                              !jyg
2984 
2985        endif ! play
2986 
2987       enddo ! l
2988
2989!       do l = 1,llm
2990!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2991!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2992!       enddo
2993 
2994          return
2995          end
2996
2997!======================================================================
2998        SUBROUTINE interp_astex_time(day,day1,annee_ref                    &
2999     &             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex         &
3000     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex        &
3001     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof   &
3002     &             ,ufa_prof,vfa_prof)
3003        implicit none
3004
3005!---------------------------------------------------------------------------------------
3006! Time interpolation of a 2D field to the timestep corresponding to day
3007!
3008! day: current julian day (e.g. 717538.2)
3009! day1: first day of the simulation
3010! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
3011! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
3012!---------------------------------------------------------------------------------------
3013
3014! inputs:
3015        integer annee_ref
3016        integer nt_astex,nlev_astex
3017        integer year_ini_astex
3018        real day, day1,day_ini_astex,dt_astex
3019        real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
3020        real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
3021! outputs:
3022        real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
3023! local:
3024        integer it_astex1, it_astex2
3025        real timeit,time_astex1,time_astex2,frac
3026
3027! Check that initial day of the simulation consistent with ASTEX period:
3028       if (annee_ref.ne.1992 ) then
3029        print*,'Pour Astex, annee_ref doit etre 1992 '
3030        print*,'Changer annee_ref dans run.def'
3031        stop
3032       endif
3033       if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
3034        print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
3035        print*,'Changer dayref dans run.def'
3036        stop
3037       endif
3038
3039! Determine timestep relative to the 1st day of TOGA-COARE:
3040!       timeit=(day-day1)*86400.
3041!       if (annee_ref.eq.1992) then
3042!        timeit=(day-day_ini_astex)*86400.
3043!       else
3044!        timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
3045!       endif
3046      timeit=(day-day_ini_astex)*86400
3047
3048! Determine the closest observation times:
3049       it_astex1=INT(timeit/dt_astex)+1
3050       it_astex2=it_astex1 + 1
3051       time_astex1=(it_astex1-1)*dt_astex
3052       time_astex2=(it_astex2-1)*dt_astex
3053       print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
3054       print *,'it_astex1,it_astex2,time_astex1,time_astex2',              &
3055     &          it_astex1,it_astex2,time_astex1,time_astex2
3056
3057       if (it_astex1 .ge. nt_astex) then
3058        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '          &
3059     &        ,day,it_astex1,it_astex2,timeit/86400.
3060        stop
3061       endif
3062
3063! time interpolation:
3064       frac=(time_astex2-timeit)/(time_astex2-time_astex1)
3065       frac=max(frac,0.0)
3066
3067       div_prof = div_astex(it_astex2)                                     &
3068     &          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
3069       ts_prof = ts_astex(it_astex2)                                        &
3070     &          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
3071       ug_prof = ug_astex(it_astex2)                                       &
3072     &          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
3073       vg_prof = vg_astex(it_astex2)                                       &
3074     &          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
3075       ufa_prof = ufa_astex(it_astex2)                                     &
3076     &          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
3077       vfa_prof = vfa_astex(it_astex2)                                     &
3078     &          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
3079
3080         print*,                                                           &
3081     &'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',       &
3082     &day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,                 &
3083     &it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
3084
3085        return
3086        END
3087
3088!======================================================================
3089        SUBROUTINE interp_toga_time(day,day1,annee_ref                     &
3090     &             ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga   &
3091     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga   &
3092     &             ,ht_toga,vt_toga,hq_toga,vq_toga                        &
3093     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof   &
3094     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
3095        implicit none
3096
3097!---------------------------------------------------------------------------------------
3098! Time interpolation of a 2D field to the timestep corresponding to day
3099!
3100! day: current julian day (e.g. 717538.2)
3101! day1: first day of the simulation
3102! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
3103! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
3104!---------------------------------------------------------------------------------------
3105
3106#include "compar1d.h"
3107
3108! inputs:
3109        integer annee_ref
3110        integer nt_toga,nlev_toga
3111        integer year_ini_toga
3112        real day, day1,day_ini_toga,dt_toga
3113        real ts_toga(nt_toga)
3114        real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
3115        real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
3116        real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
3117        real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
3118        real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
3119! outputs:
3120        real ts_prof
3121        real plev_prof(nlev_toga),t_prof(nlev_toga)
3122        real q_prof(nlev_toga),u_prof(nlev_toga)
3123        real v_prof(nlev_toga),w_prof(nlev_toga)
3124        real ht_prof(nlev_toga),vt_prof(nlev_toga)
3125        real hq_prof(nlev_toga),vq_prof(nlev_toga)
3126! local:
3127        integer it_toga1, it_toga2,k
3128        real timeit,time_toga1,time_toga2,frac
3129
3130
3131        if (forcing_type.eq.2) then
3132! Check that initial day of the simulation consistent with TOGA-COARE period:
3133       if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
3134        print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
3135        print*,'Changer annee_ref dans run.def'
3136        stop
3137       endif
3138       if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
3139        print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'
3140        print*,'Changer dayref dans run.def'
3141        stop
3142       endif
3143       if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
3144        print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
3145        print*,'Changer dayref ou nday dans run.def'
3146        stop
3147       endif
3148
3149       else if (forcing_type.eq.4) then
3150
3151! Check that initial day of the simulation consistent with TWP-ICE period:
3152       if (annee_ref.ne.2006) then
3153        print*,'Pour TWP-ICE, annee_ref doit etre 2006'
3154        print*,'Changer annee_ref dans run.def'
3155        stop
3156       endif
3157       if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
3158        print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
3159        print*,'Changer dayref dans run.def'
3160        stop
3161       endif
3162       if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
3163        print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
3164        print*,'Changer dayref ou nday dans run.def'
3165        stop
3166       endif
3167
3168       endif
3169
3170! Determine timestep relative to the 1st day of TOGA-COARE:
3171!       timeit=(day-day1)*86400.
3172!       if (annee_ref.eq.1992) then
3173!        timeit=(day-day_ini_toga)*86400.
3174!       else
3175!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3176!       endif
3177      timeit=(day-day_ini_toga)*86400
3178
3179! Determine the closest observation times:
3180       it_toga1=INT(timeit/dt_toga)+1
3181       it_toga2=it_toga1 + 1
3182       time_toga1=(it_toga1-1)*dt_toga
3183       time_toga2=(it_toga2-1)*dt_toga
3184
3185       if (it_toga1 .ge. nt_toga) then
3186        write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: '            &
3187     &        ,day,it_toga1,it_toga2,timeit/86400.
3188        stop
3189       endif
3190
3191! time interpolation:
3192       frac=(time_toga2-timeit)/(time_toga2-time_toga1)
3193       frac=max(frac,0.0)
3194
3195       ts_prof = ts_toga(it_toga2)                                         &
3196     &          -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
3197
3198!        print*,
3199!     :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
3200!     :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
3201
3202       do k=1,nlev_toga
3203        plev_prof(k) = 100.*(plev_toga(k,it_toga2)                         &
3204     &          -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
3205        t_prof(k) = t_toga(k,it_toga2)                                     &
3206     &          -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
3207        q_prof(k) = q_toga(k,it_toga2)                                     &
3208     &          -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
3209        u_prof(k) = u_toga(k,it_toga2)                                     &
3210     &          -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
3211        v_prof(k) = v_toga(k,it_toga2)                                     &
3212     &          -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
3213        w_prof(k) = w_toga(k,it_toga2)                                     &
3214     &          -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
3215        ht_prof(k) = ht_toga(k,it_toga2)                                   &
3216     &          -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
3217        vt_prof(k) = vt_toga(k,it_toga2)                                   &
3218     &          -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
3219        hq_prof(k) = hq_toga(k,it_toga2)                                   &
3220     &          -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
3221        vq_prof(k) = vq_toga(k,it_toga2)                                   &
3222     &          -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
3223        enddo
3224
3225        return
3226        END
3227
3228!======================================================================
3229        SUBROUTINE interp_dice_time(day,day1,annee_ref                    &
3230     &             ,year_ini_dice,day_ini_dice,nt_dice,dt_dice            &
3231     &             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice       &
3232     &             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice         &
3233     &             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice     &
3234     &             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof         &
3235     &             ,ustar_prof,psurf_prof,ug_prof,vg_prof                 &
3236     &             ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
3237        implicit none
3238
3239!---------------------------------------------------------------------------------------
3240! Time interpolation of a 2D field to the timestep corresponding to day
3241!
3242! day: current julian day (e.g. 717538.2)
3243! day1: first day of the simulation
3244! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
3245! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
3246!---------------------------------------------------------------------------------------
3247
3248#include "compar1d.h"
3249
3250! inputs:
3251        integer annee_ref
3252        integer nt_dice,nlev_dice
3253        integer year_ini_dice
3254        real day, day1,day_ini_dice,dt_dice
3255        real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
3256        real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
3257        real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
3258        real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
3259        real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
3260        real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
3261! outputs:
3262        real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
3263        real ustar_prof,psurf_prof,ug_prof,vg_prof
3264        real ht_prof(nlev_dice),hq_prof(nlev_dice)
3265        real hu_prof(nlev_dice),hv_prof(nlev_dice)
3266        real w_prof(nlev_dice),omega_prof(nlev_dice)
3267! local:
3268        integer it_dice1, it_dice2,k
3269        real timeit,time_dice1,time_dice2,frac
3270
3271
3272        if (forcing_type.eq.7) then
3273! Check that initial day of the simulation consistent with Dice period:
3274       print *,'annee_ref=',annee_ref
3275       print *,'day1=',day1
3276       print *,'day_ini_dice=',day_ini_dice
3277       if (annee_ref.ne.1999) then
3278        print*,'Pour Dice, annee_ref doit etre 1999'
3279        print*,'Changer annee_ref dans run.def'
3280        stop
3281       endif
3282       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
3283        print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
3284        print*,'Changer dayref dans run.def',day1,day_ini_dice
3285        stop
3286       endif
3287       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
3288        print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
3289        print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
3290        stop
3291       endif
3292
3293       endif
3294
3295! Determine timestep relative to the 1st day of TOGA-COARE:
3296!       timeit=(day-day1)*86400.
3297!       if (annee_ref.eq.1992) then
3298!        timeit=(day-day_ini_dice)*86400.
3299!       else
3300!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3301!       endif
3302      timeit=(day-day_ini_dice)*86400
3303
3304! Determine the closest observation times:
3305       it_dice1=INT(timeit/dt_dice)+1
3306       it_dice2=it_dice1 + 1
3307       time_dice1=(it_dice1-1)*dt_dice
3308       time_dice2=(it_dice2-1)*dt_dice
3309
3310       if (it_dice1 .ge. nt_dice) then
3311        write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
3312        stop
3313       endif
3314
3315! time interpolation:
3316       frac=(time_dice2-timeit)/(time_dice2-time_dice1)
3317       frac=max(frac,0.0)
3318
3319       shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
3320       lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
3321       lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
3322       swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
3323       tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
3324       ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
3325       psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
3326       ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
3327       vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
3328
3329!        print*,
3330!     :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
3331!     :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
3332
3333       do k=1,nlev_dice
3334        ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
3335        hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
3336        hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
3337        hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
3338        w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
3339        omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
3340        enddo
3341
3342        return
3343        END
3344
3345!======================================================================
3346        SUBROUTINE interp_gabls4_time(day,day1,annee_ref                              &
3347     &             ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4    &
3348     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                          &
3349     &             ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
3350        implicit none
3351
3352!---------------------------------------------------------------------------------------
3353! Time interpolation of a 2D field to the timestep corresponding to day
3354!
3355! day: current julian day
3356! day1: first day of the simulation
3357! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
3358! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
3359!---------------------------------------------------------------------------------------
3360
3361#include "compar1d.h"
3362
3363! inputs:
3364        integer annee_ref
3365        integer nt_gabls4,nlev_gabls4
3366        integer year_ini_gabls4
3367        real day, day1,day_ini_gabls4,dt_gabls4
3368        real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
3369        real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
3370        real tg_gabls4(nt_gabls4), tg_prof
3371! outputs:
3372        real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
3373        real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
3374! local:
3375        integer it_gabls41, it_gabls42,k
3376        real timeit,time_gabls41,time_gabls42,frac
3377
3378
3379
3380! Check that initial day of the simulation consistent with gabls4 period:
3381       if (forcing_type.eq.8 ) then
3382       print *,'annee_ref=',annee_ref
3383       print *,'day1=',day1
3384       print *,'day_ini_gabls4=',day_ini_gabls4
3385       if (annee_ref.ne.2009) then
3386        print*,'Pour gabls4, annee_ref doit etre 2009'
3387        print*,'Changer annee_ref dans run.def'
3388        stop
3389       endif
3390       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
3391        print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
3392        print*,'Changer dayref dans run.def',day1,day_ini_gabls4
3393        stop
3394       endif
3395       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
3396        print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
3397        print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
3398        stop
3399       endif
3400       endif
3401
3402      timeit=(day-day_ini_gabls4)*86400
3403       print *,'day,day_ini_gabls4=',day,day_ini_gabls4
3404       print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
3405
3406! Determine the closest observation times:
3407       it_gabls41=INT(timeit/dt_gabls4)+1
3408       it_gabls42=it_gabls41 + 1
3409       time_gabls41=(it_gabls41-1)*dt_gabls4
3410       time_gabls42=(it_gabls42-1)*dt_gabls4
3411
3412       if (it_gabls41 .ge. nt_gabls4) then
3413        write(*,*) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
3414        stop
3415       endif
3416
3417! time interpolation:
3418       frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
3419       frac=max(frac,0.0)
3420
3421
3422       do k=1,nlev_gabls4
3423        ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
3424        vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
3425        ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
3426        hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
3427        enddo
3428        tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
3429        return
3430        END
3431
3432!======================================================================
3433        SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
3434     &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
3435     &             ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu         &
3436     &             ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
3437        implicit none
3438
3439!---------------------------------------------------------------------------------------
3440! Time interpolation of a 2D field to the timestep corresponding to day
3441!
3442! day: current julian day (e.g. 717538.2)
3443! day1: first day of the simulation
3444! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
3445! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
3446! fs= sensible flux
3447! fl= latent flux
3448! at,rt,aqt= advective and radiative tendencies
3449!---------------------------------------------------------------------------------------
3450
3451! inputs:
3452        integer annee_ref
3453        integer nt_armcu,nlev_armcu
3454        integer year_ini_armcu
3455        real day, day1,day_ini_armcu,dt_armcu
3456        real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
3457        real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
3458! outputs:
3459        real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
3460! local:
3461        integer it_armcu1, it_armcu2,k
3462        real timeit,time_armcu1,time_armcu2,frac
3463
3464! Check that initial day of the simulation consistent with ARMCU period:
3465       if (annee_ref.ne.1997 ) then
3466        print*,'Pour ARMCU, annee_ref doit etre 1997 '
3467        print*,'Changer annee_ref dans run.def'
3468        stop
3469       endif
3470
3471      timeit=(day-day_ini_armcu)*86400
3472
3473! Determine the closest observation times:
3474       it_armcu1=INT(timeit/dt_armcu)+1
3475       it_armcu2=it_armcu1 + 1
3476       time_armcu1=(it_armcu1-1)*dt_armcu
3477       time_armcu2=(it_armcu2-1)*dt_armcu
3478       print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
3479       print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',              &
3480     &          it_armcu1,it_armcu2,time_armcu1,time_armcu2
3481
3482       if (it_armcu1 .ge. nt_armcu) then
3483        write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '          &
3484     &        ,day,it_armcu1,it_armcu2,timeit/86400.
3485        stop
3486       endif
3487
3488! time interpolation:
3489       frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
3490       frac=max(frac,0.0)
3491
3492       fs_prof = fs_armcu(it_armcu2)                                       &
3493     &          -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
3494       fl_prof = fl_armcu(it_armcu2)                                       &
3495     &          -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
3496       at_prof = at_armcu(it_armcu2)                                       &
3497     &          -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
3498       rt_prof = rt_armcu(it_armcu2)                                       &
3499     &          -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
3500       aqt_prof = aqt_armcu(it_armcu2)                                       &
3501     &          -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
3502
3503         print*,                                                           &
3504     &'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',       &
3505     &day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,                 &
3506     &it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
3507
3508        return
3509        END
3510
3511!=====================================================================
3512      subroutine readprofiles(nlev_max,kmax,ntrac,height,                  &
3513     &           thlprof,qtprof,uprof,                                     &
3514     &           vprof,e12prof,ugprof,vgprof,                              &
3515     &           wfls,dqtdxls,dqtdyls,dqtdtls,                             &
3516     &           thlpcar,tracer,nt1,nt2)
3517      implicit none
3518
3519        integer nlev_max,kmax,kmax2,ntrac
3520        logical :: llesread = .true.
3521
3522        real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),          &
3523     &       uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),            &
3524     &       ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),             &
3525     &       dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),        &
3526     &           thlpcar(nlev_max),tracer(nlev_max,ntrac)
3527
3528        real height1(nlev_max)
3529
3530        integer, parameter :: ilesfile=1
3531        integer :: ierr,k,itrac,nt1,nt2
3532
3533        if(.not.(llesread)) return
3534
3535       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3536        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3537        read (ilesfile,*) kmax
3538        do k=1,kmax
3539          read (ilesfile,*) height1(k),thlprof(k),qtprof (k),               &
3540     &                      uprof (k),vprof  (k),e12prof(k)
3541        enddo
3542        close(ilesfile)
3543
3544       open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
3545        if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
3546        read (ilesfile,*) kmax2
3547        if (kmax .ne. kmax2) then
3548          print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3549          print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3550          stop 'lecture profiles'
3551        endif
3552        do k=1,kmax
3553          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),         &
3554     &                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
3555        end do
3556        do k=1,kmax
3557          if (height(k) .ne. height1(k)) then
3558            print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3559            print *, 'les niveaux different : ',k,height1(k), height(k)
3560            stop
3561          endif
3562        end do
3563        close(ilesfile)
3564
3565       open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
3566        if (ierr /= 0) then
3567            print*,'WARNING : trac.inp does not exist'
3568        else
3569        read (ilesfile,*) kmax2,nt1,nt2
3570        if (nt2>ntrac) then
3571          stop'Augmenter le nombre de traceurs dans traceur.def'
3572        endif
3573        if (kmax .ne. kmax2) then
3574          print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3575          print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3576          stop 'lecture profiles'
3577        endif
3578        do k=1,kmax
3579          read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
3580        end do
3581        close(ilesfile)
3582        endif
3583
3584        return
3585        end
3586!======================================================================
3587      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,       &
3588     &       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
3589!======================================================================
3590      implicit none
3591
3592        integer nlev_max,kmax
3593        logical :: llesread = .true.
3594
3595        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3596        real thlprof(nlev_max)
3597        real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3598        real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
3599
3600        integer, parameter :: ilesfile=1
3601        integer :: k,ierr
3602
3603        if(.not.(llesread)) return
3604
3605       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3606        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3607        read (ilesfile,*) kmax
3608        do k=1,kmax
3609          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
3610     &                      qprof (k),uprof(k),  vprof(k),  wprof(k),      &
3611     &                      omega (k),o3mmr(k)
3612        enddo
3613        close(ilesfile)
3614
3615        return
3616        end
3617
3618!======================================================================
3619      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,       &
3620     &    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
3621!======================================================================
3622      implicit none
3623
3624        integer nlev_max,kmax
3625        logical :: llesread = .true.
3626
3627        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),             &
3628     &  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),               &
3629     &  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),                  &
3630     &  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
3631
3632        integer, parameter :: ilesfile=1
3633        integer :: ierr,k
3634
3635        if(.not.(llesread)) return
3636
3637       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3638        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3639        read (ilesfile,*) kmax
3640        do k=1,kmax
3641          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
3642     &                qvprof (k),qlprof (k),qtprof (k),                    &
3643     &                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
3644        enddo
3645        close(ilesfile)
3646
3647        return
3648        end
3649
3650
3651
3652!======================================================================
3653      subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,       &
3654     &       vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
3655!======================================================================
3656      implicit none
3657
3658        integer nlev_max,kmax
3659        logical :: llesread = .true.
3660
3661        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3662        real thetaprof(nlev_max),rvprof(nlev_max)
3663        real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3664        real aprof(nlev_max+1),bprof(nlev_max+1)
3665
3666        integer, parameter :: ilesfile=1
3667        integer, parameter :: ifile=2
3668        integer :: ierr,jtot,k
3669
3670        if(.not.(llesread)) return
3671
3672! Read profiles at full levels
3673       IF(nlev_max.EQ.19) THEN
3674       open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
3675       print *,'On ouvre prof.inp.19'
3676       ELSE
3677       open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
3678       print *,'On ouvre prof.inp.40'
3679       ENDIF
3680        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3681        read (ilesfile,*) kmax
3682        do k=1,kmax
3683          read (ilesfile,*) height(k)    ,pprof(k),  uprof(k), vprof(k),   &
3684     &                      thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
3685        enddo
3686        close(ilesfile)
3687
3688! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf) 
3689       IF(nlev_max.EQ.19) THEN
3690       open (ifile,file='proh.inp.19',status='old',iostat=ierr)
3691       print *,'On ouvre proh.inp.19'
3692       if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
3693       ELSE
3694       open (ifile,file='proh.inp.40',status='old',iostat=ierr)
3695       print *,'On ouvre proh.inp.40'
3696       if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
3697       ENDIF
3698        read (ifile,*) kmax
3699        do k=1,kmax
3700          read (ifile,*) jtot,aprof(k),bprof(k)
3701        enddo
3702        close(ifile)
3703
3704        return
3705        end
3706
3707!=====================================================================
3708      subroutine read_fire(fich_fire,nlevel,ntime                          &
3709     &     ,zz,thl,qt,u,v,tke                                              &
3710     &     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3711
3712!program reading forcings of the FIRE case study
3713
3714
3715      implicit none
3716
3717#include "netcdf.inc"
3718
3719      integer ntime,nlevel
3720      character*80 :: fich_fire
3721      real*8 zz(nlevel)
3722
3723      real*8 thl(nlevel)
3724      real*8 qt(nlevel),u(nlevel)
3725      real*8 v(nlevel),tke(nlevel)
3726      real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3727      real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3728      real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3729
3730      integer nid, ierr
3731      integer nbvar3d
3732      parameter(nbvar3d=30)
3733      integer var3didin(nbvar3d)
3734
3735      ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3736      if (ierr.NE.NF_NOERR) then
3737         write(*,*) 'ERROR: Pb opening forcings nc file '
3738         write(*,*) NF_STRERROR(ierr)
3739         stop ""
3740      endif
3741
3742
3743       ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) 
3744         if(ierr/=NF_NOERR) then
3745           write(*,*) NF_STRERROR(ierr)
3746           stop 'lev'
3747         endif
3748
3749
3750      ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3751         if(ierr/=NF_NOERR) then
3752           write(*,*) NF_STRERROR(ierr)
3753           stop 'temp'
3754         endif
3755
3756      ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3757         if(ierr/=NF_NOERR) then
3758           write(*,*) NF_STRERROR(ierr)
3759           stop 'qv'
3760         endif
3761
3762      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3763         if(ierr/=NF_NOERR) then
3764           write(*,*) NF_STRERROR(ierr)
3765           stop 'u'
3766         endif
3767
3768      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3769         if(ierr/=NF_NOERR) then
3770           write(*,*) NF_STRERROR(ierr)
3771           stop 'v'
3772         endif
3773
3774      ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3775         if(ierr/=NF_NOERR) then
3776           write(*,*) NF_STRERROR(ierr)
3777           stop 'tke'
3778         endif
3779
3780      ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3781         if(ierr/=NF_NOERR) then
3782           write(*,*) NF_STRERROR(ierr)
3783           stop 'ug'
3784         endif
3785
3786      ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3787         if(ierr/=NF_NOERR) then
3788           write(*,*) NF_STRERROR(ierr)
3789           stop 'vg'
3790         endif
3791     
3792      ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3793         if(ierr/=NF_NOERR) then
3794           write(*,*) NF_STRERROR(ierr)
3795           stop 'wls'
3796         endif
3797
3798      ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3799         if(ierr/=NF_NOERR) then
3800           write(*,*) NF_STRERROR(ierr)
3801           stop 'dqtdx'
3802         endif
3803
3804      ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3805         if(ierr/=NF_NOERR) then
3806           write(*,*) NF_STRERROR(ierr)
3807           stop 'dqtdy'
3808      endif
3809
3810      ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3811         if(ierr/=NF_NOERR) then
3812           write(*,*) NF_STRERROR(ierr)
3813           stop 'dqtdt'
3814      endif
3815
3816      ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3817         if(ierr/=NF_NOERR) then
3818           write(*,*) NF_STRERROR(ierr)
3819           stop 'thl_rad'
3820      endif
3821!dimensions lecture
3822!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
3823 
3824#ifdef NC_DOUBLE
3825         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3826#else
3827         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3828#endif
3829         if(ierr/=NF_NOERR) then
3830            write(*,*) NF_STRERROR(ierr)
3831            stop "getvarup"
3832         endif
3833!          write(*,*)'lecture z ok',zz
3834
3835#ifdef NC_DOUBLE
3836         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3837#else
3838         ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3839#endif
3840         if(ierr/=NF_NOERR) then
3841            write(*,*) NF_STRERROR(ierr)
3842            stop "getvarup"
3843         endif
3844!          write(*,*)'lecture thl ok',thl
3845
3846#ifdef NC_DOUBLE
3847         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3848#else
3849         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3850#endif
3851         if(ierr/=NF_NOERR) then
3852            write(*,*) NF_STRERROR(ierr)
3853            stop "getvarup"
3854         endif
3855!          write(*,*)'lecture qt ok',qt
3856 
3857#ifdef NC_DOUBLE
3858         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3859#else
3860         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3861#endif
3862         if(ierr/=NF_NOERR) then
3863            write(*,*) NF_STRERROR(ierr)
3864            stop "getvarup"
3865         endif
3866!          write(*,*)'lecture u ok',u
3867
3868#ifdef NC_DOUBLE
3869         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3870#else
3871         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3872#endif
3873         if(ierr/=NF_NOERR) then
3874            write(*,*) NF_STRERROR(ierr)
3875            stop "getvarup"
3876         endif
3877!          write(*,*)'lecture v ok',v
3878
3879#ifdef NC_DOUBLE
3880         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3881#else
3882         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3883#endif
3884         if(ierr/=NF_NOERR) then
3885            write(*,*) NF_STRERROR(ierr)
3886            stop "getvarup"
3887         endif
3888!          write(*,*)'lecture tke ok',tke
3889
3890#ifdef NC_DOUBLE
3891         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3892#else
3893         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3894#endif
3895         if(ierr/=NF_NOERR) then
3896            write(*,*) NF_STRERROR(ierr)
3897            stop "getvarup"
3898         endif
3899!          write(*,*)'lecture ug ok',ug
3900
3901#ifdef NC_DOUBLE
3902         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3903#else
3904         ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3905#endif
3906         if(ierr/=NF_NOERR) then
3907            write(*,*) NF_STRERROR(ierr)
3908            stop "getvarup"
3909         endif
3910!          write(*,*)'lecture vg ok',vg
3911
3912#ifdef NC_DOUBLE
3913         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3914#else
3915         ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3916#endif
3917         if(ierr/=NF_NOERR) then
3918            write(*,*) NF_STRERROR(ierr)
3919            stop "getvarup"
3920         endif
3921!          write(*,*)'lecture wls ok',wls
3922
3923#ifdef NC_DOUBLE
3924         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3925#else
3926         ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3927#endif
3928         if(ierr/=NF_NOERR) then
3929            write(*,*) NF_STRERROR(ierr)
3930            stop "getvarup"
3931         endif
3932!          write(*,*)'lecture dqtdx ok',dqtdx
3933
3934#ifdef NC_DOUBLE
3935         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3936#else
3937         ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3938#endif
3939         if(ierr/=NF_NOERR) then
3940            write(*,*) NF_STRERROR(ierr)
3941            stop "getvarup"
3942         endif
3943!          write(*,*)'lecture dqtdy ok',dqtdy
3944
3945#ifdef NC_DOUBLE
3946         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3947#else
3948         ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3949#endif
3950         if(ierr/=NF_NOERR) then
3951            write(*,*) NF_STRERROR(ierr)
3952            stop "getvarup"
3953         endif
3954!          write(*,*)'lecture dqtdt ok',dqtdt
3955
3956#ifdef NC_DOUBLE
3957         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3958#else
3959         ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3960#endif
3961         if(ierr/=NF_NOERR) then
3962            write(*,*) NF_STRERROR(ierr)
3963            stop "getvarup"
3964         endif
3965!          write(*,*)'lecture thl_rad ok',thl_rad
3966
3967         return 
3968         end subroutine read_fire
3969!=====================================================================
3970      subroutine read_dice(fich_dice,nlevel,ntime                         &
3971     &     ,zz,pres,t,qv,u,v,o3                                          &
3972     &     ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg                        &
3973     &     ,hadvt,hadvq,hadvu,hadvv,w,omega)
3974
3975!program reading initial profils and forcings of the Dice case study
3976
3977
3978      implicit none
3979
3980#include "netcdf.inc"
3981#include "YOMCST.h"
3982
3983      integer ntime,nlevel
3984      integer l,k
3985      character*80 :: fich_dice
3986      real*8 time(ntime)
3987      real*8 zz(nlevel)
3988
3989      real*8 th(nlevel),pres(nlevel),t(nlevel)
3990      real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
3991      real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
3992      real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
3993      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
3994      real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
3995      real*8 pzero
3996
3997      integer nid, ierr
3998      integer nbvar3d
3999      parameter(nbvar3d=30)
4000      integer var3didin(nbvar3d)
4001
4002      pzero=100000.
4003      ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
4004      if (ierr.NE.NF_NOERR) then
4005         write(*,*) 'ERROR: Pb opening forcings nc file '
4006         write(*,*) NF_STRERROR(ierr)
4007         stop ""
4008      endif
4009
4010
4011       ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 
4012         if(ierr/=NF_NOERR) then
4013           write(*,*) NF_STRERROR(ierr)
4014           stop 'height'
4015         endif
4016
4017       ierr=NF_INQ_VARID(nid,"pf",var3didin(11)) 
4018         if(ierr/=NF_NOERR) then
4019           write(*,*) NF_STRERROR(ierr)
4020           stop 'pf'
4021         endif
4022
4023      ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
4024         if(ierr/=NF_NOERR) then
4025           write(*,*) NF_STRERROR(ierr)
4026           stop 'theta'
4027         endif
4028
4029      ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
4030         if(ierr/=NF_NOERR) then
4031           write(*,*) NF_STRERROR(ierr)
4032           stop 'qv'
4033         endif
4034
4035      ierr=NF_INQ_VARID(nid,"u",var3didin(14))
4036         if(ierr/=NF_NOERR) then
4037           write(*,*) NF_STRERROR(ierr)
4038           stop 'u'
4039         endif
4040
4041      ierr=NF_INQ_VARID(nid,"v",var3didin(15))
4042         if(ierr/=NF_NOERR) then
4043           write(*,*) NF_STRERROR(ierr)
4044           stop 'v'
4045         endif
4046
4047      ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
4048         if(ierr/=NF_NOERR) then
4049           write(*,*) NF_STRERROR(ierr)
4050           stop 'o3'
4051         endif
4052
4053      ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
4054         if(ierr/=NF_NOERR) then
4055           write(*,*) NF_STRERROR(ierr)
4056           stop 'shf'
4057         endif
4058
4059      ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
4060         if(ierr/=NF_NOERR) then
4061           write(*,*) NF_STRERROR(ierr)
4062           stop 'lhf'
4063         endif
4064     
4065      ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
4066         if(ierr/=NF_NOERR) then
4067           write(*,*) NF_STRERROR(ierr)
4068           stop 'lwup'
4069         endif
4070
4071      ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
4072         if(ierr/=NF_NOERR) then
4073           write(*,*) NF_STRERROR(ierr)
4074           stop 'dqtdx'
4075         endif
4076
4077      ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
4078         if(ierr/=NF_NOERR) then
4079           write(*,*) NF_STRERROR(ierr)
4080           stop 'Tg'
4081      endif
4082
4083      ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
4084         if(ierr/=NF_NOERR) then
4085           write(*,*) NF_STRERROR(ierr)
4086           stop 'ustar'
4087      endif
4088
4089      ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
4090         if(ierr/=NF_NOERR) then
4091           write(*,*) NF_STRERROR(ierr)
4092           stop 'psurf'
4093      endif
4094
4095      ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
4096         if(ierr/=NF_NOERR) then
4097           write(*,*) NF_STRERROR(ierr)
4098           stop 'Ug'
4099      endif
4100
4101      ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
4102         if(ierr/=NF_NOERR) then
4103           write(*,*) NF_STRERROR(ierr)
4104           stop 'Vg'
4105      endif
4106
4107      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
4108         if(ierr/=NF_NOERR) then
4109           write(*,*) NF_STRERROR(ierr)
4110           stop 'hadvT'
4111      endif
4112
4113      ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
4114         if(ierr/=NF_NOERR) then
4115           write(*,*) NF_STRERROR(ierr)
4116           stop 'hadvq'
4117      endif
4118
4119      ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
4120         if(ierr/=NF_NOERR) then
4121           write(*,*) NF_STRERROR(ierr)
4122           stop 'hadvu'
4123      endif
4124
4125      ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
4126         if(ierr/=NF_NOERR) then
4127           write(*,*) NF_STRERROR(ierr)
4128           stop 'hadvv'
4129      endif
4130
4131      ierr=NF_INQ_VARID(nid,"w",var3didin(21))
4132         if(ierr/=NF_NOERR) then
4133           write(*,*) NF_STRERROR(ierr)
4134           stop 'w'
4135      endif
4136
4137      ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
4138         if(ierr/=NF_NOERR) then
4139           write(*,*) NF_STRERROR(ierr)
4140           stop 'omega'
4141      endif
4142!dimensions lecture
4143!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
4144 
4145#ifdef NC_DOUBLE
4146         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
4147#else
4148         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
4149#endif
4150         if(ierr/=NF_NOERR) then
4151            write(*,*) NF_STRERROR(ierr)
4152            stop "getvarup"
4153         endif
4154!          write(*,*)'lecture zz ok',zz
4155 
4156#ifdef NC_DOUBLE
4157         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
4158#else
4159         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
4160#endif
4161         if(ierr/=NF_NOERR) then
4162            write(*,*) NF_STRERROR(ierr)
4163            stop "getvarup"
4164         endif
4165!          write(*,*)'lecture pres ok',pres
4166
4167#ifdef NC_DOUBLE
4168         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
4169#else
4170         ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
4171#endif
4172         if(ierr/=NF_NOERR) then
4173            write(*,*) NF_STRERROR(ierr)
4174            stop "getvarup"
4175         endif
4176!          write(*,*)'lecture th ok',th
4177           do k=1,nlevel
4178             t(k)=th(k)*(pres(k)/pzero)**rkappa
4179           enddo
4180
4181#ifdef NC_DOUBLE
4182         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
4183#else
4184         ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
4185#endif
4186         if(ierr/=NF_NOERR) then
4187            write(*,*) NF_STRERROR(ierr)
4188            stop "getvarup"
4189         endif
4190!          write(*,*)'lecture qv ok',qv
4191 
4192#ifdef NC_DOUBLE
4193         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
4194#else
4195         ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
4196#endif
4197         if(ierr/=NF_NOERR) then
4198            write(*,*) NF_STRERROR(ierr)
4199            stop "getvarup"
4200         endif
4201!          write(*,*)'lecture u ok',u
4202
4203#ifdef NC_DOUBLE
4204         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
4205#else
4206         ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
4207#endif
4208         if(ierr/=NF_NOERR) then
4209            write(*,*) NF_STRERROR(ierr)
4210            stop "getvarup"
4211         endif
4212!          write(*,*)'lecture v ok',v
4213
4214#ifdef NC_DOUBLE
4215         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
4216#else
4217         ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
4218#endif
4219         if(ierr/=NF_NOERR) then
4220            write(*,*) NF_STRERROR(ierr)
4221            stop "getvarup"
4222         endif
4223!          write(*,*)'lecture o3 ok',o3
4224
4225#ifdef NC_DOUBLE
4226         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
4227#else
4228         ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
4229#endif
4230         if(ierr/=NF_NOERR) then
4231            write(*,*) NF_STRERROR(ierr)
4232            stop "getvarup"
4233         endif
4234!          write(*,*)'lecture shf ok',shf
4235
4236#ifdef NC_DOUBLE
4237         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
4238#else
4239         ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
4240#endif
4241         if(ierr/=NF_NOERR) then
4242            write(*,*) NF_STRERROR(ierr)
4243            stop "getvarup"
4244         endif
4245!          write(*,*)'lecture lhf ok',lhf
4246
4247#ifdef NC_DOUBLE
4248         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
4249#else
4250         ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
4251#endif
4252         if(ierr/=NF_NOERR) then
4253            write(*,*) NF_STRERROR(ierr)
4254            stop "getvarup"
4255         endif
4256!          write(*,*)'lecture lwup ok',lwup
4257
4258#ifdef NC_DOUBLE
4259         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
4260#else
4261         ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
4262#endif
4263         if(ierr/=NF_NOERR) then
4264            write(*,*) NF_STRERROR(ierr)
4265            stop "getvarup"
4266         endif
4267!          write(*,*)'lecture swup ok',swup
4268
4269#ifdef NC_DOUBLE
4270         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
4271#else
4272         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
4273#endif
4274         if(ierr/=NF_NOERR) then
4275            write(*,*) NF_STRERROR(ierr)
4276            stop "getvarup"
4277         endif
4278!          write(*,*)'lecture tg ok',tg
4279
4280#ifdef NC_DOUBLE
4281         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
4282#else
4283         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
4284#endif
4285         if(ierr/=NF_NOERR) then
4286            write(*,*) NF_STRERROR(ierr)
4287            stop "getvarup"
4288         endif
4289!          write(*,*)'lecture ustar ok',ustar
4290
4291#ifdef NC_DOUBLE
4292         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
4293#else
4294         ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
4295#endif
4296         if(ierr/=NF_NOERR) then
4297            write(*,*) NF_STRERROR(ierr)
4298            stop "getvarup"
4299         endif
4300!          write(*,*)'lecture psurf ok',psurf
4301
4302#ifdef NC_DOUBLE
4303         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
4304#else
4305         ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
4306#endif
4307         if(ierr/=NF_NOERR) then
4308            write(*,*) NF_STRERROR(ierr)
4309            stop "getvarup"
4310         endif
4311!          write(*,*)'lecture ug ok',ug
4312
4313#ifdef NC_DOUBLE
4314         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
4315#else
4316         ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
4317#endif
4318         if(ierr/=NF_NOERR) then
4319            write(*,*) NF_STRERROR(ierr)
4320            stop "getvarup"
4321         endif
4322!          write(*,*)'lecture vg ok',vg
4323
4324#ifdef NC_DOUBLE
4325         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
4326#else
4327         ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
4328#endif
4329         if(ierr/=NF_NOERR) then
4330            write(*,*) NF_STRERROR(ierr)
4331            stop "getvarup"
4332         endif
4333!          write(*,*)'lecture hadvt ok',hadvt
4334
4335#ifdef NC_DOUBLE
4336         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
4337#else
4338         ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
4339#endif
4340         if(ierr/=NF_NOERR) then
4341            write(*,*) NF_STRERROR(ierr)
4342            stop "getvarup"
4343         endif
4344!          write(*,*)'lecture hadvq ok',hadvq
4345
4346#ifdef NC_DOUBLE
4347         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
4348#else
4349         ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
4350#endif
4351         if(ierr/=NF_NOERR) then
4352            write(*,*) NF_STRERROR(ierr)
4353            stop "getvarup"
4354         endif
4355!          write(*,*)'lecture hadvu ok',hadvu
4356
4357#ifdef NC_DOUBLE
4358         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
4359#else
4360         ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
4361#endif
4362         if(ierr/=NF_NOERR) then
4363            write(*,*) NF_STRERROR(ierr)
4364            stop "getvarup"
4365         endif
4366!          write(*,*)'lecture hadvv ok',hadvv
4367
4368#ifdef NC_DOUBLE
4369         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
4370#else
4371         ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
4372#endif
4373         if(ierr/=NF_NOERR) then
4374            write(*,*) NF_STRERROR(ierr)
4375            stop "getvarup"
4376         endif
4377!          write(*,*)'lecture w ok',w
4378
4379#ifdef NC_DOUBLE
4380         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
4381#else
4382         ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
4383#endif
4384         if(ierr/=NF_NOERR) then
4385            write(*,*) NF_STRERROR(ierr)
4386            stop "getvarup"
4387         endif
4388!          write(*,*)'lecture omega ok',omega
4389
4390         return 
4391         end subroutine read_dice
4392!=====================================================================
4393      subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol                    &
4394     &     ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
4395
4396!program reading initial profils and forcings of the Gabls4 case study
4397
4398
4399      implicit none
4400
4401#include "netcdf.inc"
4402
4403      integer ntime,nlevel,nsol
4404      integer l,k
4405      character*80 :: fich_gabls4
4406      real*8 time(ntime)
4407
4408!  ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees
4409! dans un ordre inverse par rapport a la convention LMDZ
4410! ==> il faut tout inverser  (MPL 20141024)
4411! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
4412      real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
4413      real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
4414      real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
4415
4416      real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
4417      real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
4418      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
4419
4420      real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
4421      real*8 tg(ntime)
4422      integer nid, ierr
4423      integer nbvar3d
4424      parameter(nbvar3d=30)
4425      integer var3didin(nbvar3d)
4426
4427      ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
4428      if (ierr.NE.NF_NOERR) then
4429         write(*,*) 'ERROR: Pb opening forcings nc file '
4430         write(*,*) NF_STRERROR(ierr)
4431         stop ""
4432      endif
4433
4434
4435       ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 
4436         if(ierr/=NF_NOERR) then
4437           write(*,*) NF_STRERROR(ierr)
4438           stop 'height'
4439         endif
4440
4441      ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
4442         if(ierr/=NF_NOERR) then
4443           write(*,*) NF_STRERROR(ierr)
4444           stop 'depth_sn'
4445      endif
4446
4447      ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
4448         if(ierr/=NF_NOERR) then
4449           write(*,*) NF_STRERROR(ierr)
4450           stop 'Ug'
4451      endif
4452
4453      ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
4454         if(ierr/=NF_NOERR) then
4455           write(*,*) NF_STRERROR(ierr)
4456           stop 'Vg'
4457      endif
4458       ierr=NF_INQ_VARID(nid,"pf",var3didin(5)) 
4459         if(ierr/=NF_NOERR) then
4460           write(*,*) NF_STRERROR(ierr)
4461           stop 'pf'
4462         endif
4463
4464      ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
4465         if(ierr/=NF_NOERR) then
4466           write(*,*) NF_STRERROR(ierr)
4467           stop 'theta'
4468         endif
4469
4470      ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
4471         if(ierr/=NF_NOERR) then
4472           write(*,*) NF_STRERROR(ierr)
4473           stop 'tempe'
4474         endif
4475
4476      ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
4477         if(ierr/=NF_NOERR) then
4478           write(*,*) NF_STRERROR(ierr)
4479           stop 'qv'
4480         endif
4481
4482      ierr=NF_INQ_VARID(nid,"u",var3didin(9))
4483         if(ierr/=NF_NOERR) then
4484           write(*,*) NF_STRERROR(ierr)
4485           stop 'u'
4486         endif
4487
4488      ierr=NF_INQ_VARID(nid,"v",var3didin(10))
4489         if(ierr/=NF_NOERR) then
4490           write(*,*) NF_STRERROR(ierr)
4491           stop 'v'
4492         endif
4493
4494      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
4495         if(ierr/=NF_NOERR) then
4496           write(*,*) NF_STRERROR(ierr)
4497           stop 'hadvt'
4498         endif
4499
4500      ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
4501         if(ierr/=NF_NOERR) then
4502           write(*,*) NF_STRERROR(ierr)
4503           stop 'hadvq'
4504      endif
4505
4506      ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
4507         if(ierr/=NF_NOERR) then
4508           write(*,*) NF_STRERROR(ierr)
4509           stop 'tsnow'
4510      endif
4511
4512      ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
4513         if(ierr/=NF_NOERR) then
4514           write(*,*) NF_STRERROR(ierr)
4515           stop 'snow_density'
4516      endif
4517
4518      ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
4519         if(ierr/=NF_NOERR) then
4520           write(*,*) NF_STRERROR(ierr)
4521           stop 'Tg'
4522      endif
4523
4524
4525!dimensions lecture
4526!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
4527 
4528#ifdef NC_DOUBLE
4529         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)
4530#else
4531         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)
4532#endif
4533         if(ierr/=NF_NOERR) then
4534            write(*,*) NF_STRERROR(ierr)
4535            stop "getvarup"
4536         endif
4537 
4538#ifdef NC_DOUBLE
4539         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)
4540#else
4541         ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)
4542#endif
4543         if(ierr/=NF_NOERR) then
4544            write(*,*) NF_STRERROR(ierr)
4545            stop "getvarup"
4546         endif
4547 
4548#ifdef NC_DOUBLE
4549         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)
4550#else
4551         ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)
4552#endif
4553         if(ierr/=NF_NOERR) then
4554            write(*,*) NF_STRERROR(ierr)
4555            stop "getvarup"
4556         endif
4557 
4558#ifdef NC_DOUBLE
4559         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)
4560#else
4561         ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)
4562#endif
4563         if(ierr/=NF_NOERR) then
4564            write(*,*) NF_STRERROR(ierr)
4565            stop "getvarup"
4566         endif
4567 
4568#ifdef NC_DOUBLE
4569         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)
4570#else
4571         ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)
4572#endif
4573         if(ierr/=NF_NOERR) then
4574            write(*,*) NF_STRERROR(ierr)
4575            stop "getvarup"
4576         endif
4577
4578#ifdef NC_DOUBLE
4579         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)
4580#else
4581         ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)
4582#endif
4583         if(ierr/=NF_NOERR) then
4584            write(*,*) NF_STRERROR(ierr)
4585            stop "getvarup"
4586         endif
4587
4588#ifdef NC_DOUBLE
4589         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)
4590#else
4591         ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)
4592#endif
4593         if(ierr/=NF_NOERR) then
4594            write(*,*) NF_STRERROR(ierr)
4595            stop "getvarup"
4596         endif
4597
4598#ifdef NC_DOUBLE
4599         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)
4600#else
4601         ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)
4602#endif
4603         if(ierr/=NF_NOERR) then
4604            write(*,*) NF_STRERROR(ierr)
4605            stop "getvarup"
4606         endif
4607 
4608#ifdef NC_DOUBLE
4609         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)
4610#else
4611         ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)
4612#endif
4613         if(ierr/=NF_NOERR) then
4614            write(*,*) NF_STRERROR(ierr)
4615            stop "getvarup"
4616         endif
4617 
4618#ifdef NC_DOUBLE
4619         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)
4620#else
4621         ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)
4622#endif
4623         if(ierr/=NF_NOERR) then
4624            write(*,*) NF_STRERROR(ierr)
4625            stop "getvarup"
4626         endif
4627 
4628#ifdef NC_DOUBLE
4629         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)
4630#else
4631         ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)
4632#endif
4633         if(ierr/=NF_NOERR) then
4634            write(*,*) NF_STRERROR(ierr)
4635            stop "getvarup"
4636         endif
4637 
4638#ifdef NC_DOUBLE
4639         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)
4640#else
4641         ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)
4642#endif
4643         if(ierr/=NF_NOERR) then
4644            write(*,*) NF_STRERROR(ierr)
4645            stop "getvarup"
4646         endif
4647 
4648#ifdef NC_DOUBLE
4649         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)
4650#else
4651         ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)
4652#endif
4653         if(ierr/=NF_NOERR) then
4654            write(*,*) NF_STRERROR(ierr)
4655            stop "getvarup"
4656         endif
4657 
4658#ifdef NC_DOUBLE
4659         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)
4660#else
4661         ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)
4662#endif
4663         if(ierr/=NF_NOERR) then
4664            write(*,*) NF_STRERROR(ierr)
4665            stop "getvarup"
4666         endif
4667
4668#ifdef NC_DOUBLE
4669         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)
4670#else
4671         ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)
4672#endif
4673         if(ierr/=NF_NOERR) then
4674            write(*,*) NF_STRERROR(ierr)
4675            stop "getvarup"
4676         endif
4677
4678! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
4679         do k=1,nlevel
4680           zz(k)=zz_i(nlevel+1-k)
4681           ug(k,:)=ug_i(nlevel+1-k,:)
4682           vg(k,:)=vg_i(nlevel+1-k,:)
4683           pf(k)=pf_i(nlevel+1-k)
4684           print *,'pf=',pf(k)
4685           th(k)=th_i(nlevel+1-k)
4686           t(k)=t_i(nlevel+1-k)
4687           qv(k)=qv_i(nlevel+1-k)
4688           u(k)=u_i(nlevel+1-k)
4689           v(k)=v_i(nlevel+1-k)
4690           hadvt(k,:)=hadvt_i(nlevel+1-k,:)
4691           hadvq(k,:)=hadvq_i(nlevel+1-k,:)
4692         enddo
4693         return 
4694 end subroutine read_gabls4
4695!=====================================================================
4696
4697!     Reads CIRC input files     
4698
4699      SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
4700     
4701      parameter (ncm_1=49180)
4702#include "YOMCST.h"
4703
4704      real albsfc(ncm_1), albsfc_w(ncm_1)
4705      real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
4706           reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
4707      real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
4708      real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
4709      real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
4710      real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
4711           o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
4712!     za= zenital angle
4713!     sza= cosinus angle zenital
4714      real wavn(ncm_1), ssf(ncm_1),za,sza
4715      integer nlev
4716
4717
4718!     Open the files
4719
4720      open (11, file='Tsfc_sza_nlev_case.txt', status='old')
4721      open (12, file='level_input_case.txt', status='old')
4722      open (13, file='layer_input_case.txt', status='old')
4723      open (14, file='aerosol_input_case.txt', status='old')
4724      open (15, file='cloud_input_case.txt', status='old')
4725      open (16, file='sfcalbedo_input_case.txt', status='old')
4726     
4727!     Read scalar information
4728      do iskip=1,5
4729         read (11, *)
4730      enddo
4731      read (11, '(i8)') nlev
4732      read (11, '(f10.2)') tsfc
4733      read (11, '(f10.2)') za
4734      read (11, '(f10.4)') sw_dn_toa
4735      sza=cos(za/180.*RPI)
4736      print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
4737      close(11)
4738
4739!     Read level information
4740      read (12, *)
4741      do il=1,nlev
4742         read (12, 302) ilev, z(il), p(il), t(il)
4743         z(il)=z(il)*1000.    ! z donne en km
4744         p(il)=p(il)*100.     ! p donne en mb
4745      enddo
4746302   format (i8, f8.3, 2f9.2)
4747      close(12)
4748
4749!     Read layer information (midpoint values)
4750      do iskip=1,3
4751         read (13, *)
4752      enddo
4753      do il=1,nlev-1
4754         read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
4755                        n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
4756                        f11(il),f12(il)
4757         pm(il)=pm(il)*100.
4758      enddo
4759303   format (i8, 2f9.2, 10(2x,e13.7))     
4760      close(13)
4761     
4762!     Read aerosol layer information
4763      do iskip=1,3
4764         read (14, *)
4765      enddo
4766      read (14, '(f10.2)') aer_alpha
4767      read (14, *)
4768      read (14, *)
4769      do il=1,nlev-1
4770         read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
4771      enddo
4772304   format (i8, f9.5, 2f8.3)
4773      close(14)
4774     
4775!     Read cloud information
4776      do iskip=1,3
4777         read (15, *)
4778      enddo
4779      do il=1,nlev-1
4780         read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
4781         lwp(il)=lwp(il)/1000.          ! lwp donne en g/kg
4782         iwp(il)=iwp(il)/1000.          ! iwp donne en g/kg
4783         reliq(il)=reliq(il)/1000000.   ! reliq donne en microns
4784         reice(il)=reice(il)/1000000.   ! reice donne en microns
4785      enddo
4786305   format (i8, f8.3, 4f9.2)
4787      close(15)
4788
4789!     Read surface albedo (weighted & unweighted) and spectral solar irradiance
4790      do iskip=1,6
4791         read (16, *)
4792      enddo
4793      do icm_1=1,ncm_1
4794         read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
4795      enddo
4796306   format(f10.1, 2f12.5, f14.8)
4797      close(16)
4798 
4799      return 
4800      end subroutine read_circ
4801!=====================================================================
4802!     Reads RTMIP input files     
4803
4804      SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
4805     
4806#include "YOMCST.h"
4807
4808      real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
4809      real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
4810      integer nlev
4811
4812
4813!     Open the files
4814
4815      open (11, file='low_resolution_profile.txt', status='old')
4816     
4817!     Read level information
4818      read (11, *)
4819      do il=1,nlev_rtmip
4820         read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
4821      enddo
4822      do il=1,nlev_rtmip
4823         play(il)=pt(nlev_rtmip-il+1)*100.     ! p donne en mb
4824         temp(il)=t(nlev_rtmip-il+1)
4825         ovap(il)=h2o(nlev_rtmip-il+1)
4826         oz(il)=o3(nlev_rtmip-il+1)
4827      enddo
4828      do il=1,39
4829         plev(il)=play(il)+(play(il+1)-play(il))/2.
4830         print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
4831      enddo
4832      plev(41)=101300.
4833302   format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
4834      close(12)
4835 
4836      return 
4837      end subroutine read_rtmip
4838!=====================================================================
4839
4840!  Subroutines for nudging
4841
4842      Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
4843! ========================================================
4844  USE dimphy
4845
4846  implicit none
4847
4848! ========================================================
4849      REAL paprs(klon,klevp1)
4850      REAL pplay(klon,klev)
4851!
4852!      Variables d'etat
4853      REAL t(klon,klev)
4854      REAL q(klon,klev)
4855!
4856!   Profiles cible
4857      REAL t_targ(klon,klev)
4858      REAL rh_targ(klon,klev)
4859!
4860   INTEGER k,i
4861   REAL zx_qs
4862
4863! Declaration des constantes et des fonctions thermodynamiques
4864!
4865include "YOMCST.h"
4866include "YOETHF.h"
4867!
4868!  ----------------------------------------
4869!  Statement functions
4870include "FCTTRE.h"
4871!  ----------------------------------------
4872!
4873        DO k = 1,klev
4874         DO i = 1,klon
4875           t_targ(i,k) = t(i,k)
4876           IF (t(i,k).LT.RTT) THEN
4877              zx_qs = qsats(t(i,k))/(pplay(i,k))
4878           ELSE
4879              zx_qs = qsatl(t(i,k))/(pplay(i,k))
4880           ENDIF
4881           rh_targ(i,k) = q(i,k)/zx_qs
4882         ENDDO
4883        ENDDO
4884      print *, 't_targ',t_targ
4885      print *, 'rh_targ',rh_targ
4886!
4887!
4888      RETURN
4889      END
4890
4891      Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
4892! ========================================================
4893  USE dimphy
4894
4895  implicit none
4896
4897! ========================================================
4898      REAL paprs(klon,klevp1)
4899      REAL pplay(klon,klev)
4900!
4901!      Variables d'etat
4902      REAL u(klon,klev)
4903      REAL v(klon,klev)
4904!
4905!   Profiles cible
4906      REAL u_targ(klon,klev)
4907      REAL v_targ(klon,klev)
4908!
4909   INTEGER k,i
4910!
4911        DO k = 1,klev
4912         DO i = 1,klon
4913           u_targ(i,k) = u(i,k)
4914           v_targ(i,k) = v(i,k)
4915         ENDDO
4916        ENDDO
4917      print *, 'u_targ',u_targ
4918      print *, 'v_targ',v_targ
4919!
4920!
4921      RETURN
4922      END
4923
4924      Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q,          &
4925     &                      d_t,d_q)
4926! ========================================================
4927  USE dimphy
4928
4929  implicit none
4930
4931! ========================================================
4932      REAL dtime
4933      REAL paprs(klon,klevp1)
4934      REAL pplay(klon,klev)
4935!
4936!      Variables d'etat
4937      REAL t(klon,klev)
4938      REAL q(klon,klev)
4939!
4940! Tendances
4941      REAL d_t(klon,klev)
4942      REAL d_q(klon,klev)
4943!
4944!   Profiles cible
4945      REAL t_targ(klon,klev)
4946      REAL rh_targ(klon,klev)
4947!
4948!   Temps de relaxation
4949      REAL tau
4950!c      DATA tau /3600./
4951!!      DATA tau /5400./
4952      DATA tau /1800./
4953!
4954   INTEGER k,i
4955   REAL zx_qs, rh, tnew, d_rh, rhnew
4956
4957! Declaration des constantes et des fonctions thermodynamiques
4958!
4959include "YOMCST.h"
4960include "YOETHF.h"
4961!
4962!  ----------------------------------------
4963!  Statement functions
4964include "FCTTRE.h"
4965!  ----------------------------------------
4966!
4967        print *,'dtime, tau ',dtime,tau
4968        print *, 't_targ',t_targ
4969        print *, 'rh_targ',rh_targ
4970        print *,'temp ',t
4971        print *,'hum ',q
4972!
4973        DO k = 1,klev
4974         DO i = 1,klon
4975           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
4976            IF (t(i,k).LT.RTT) THEN
4977               zx_qs = qsats(t(i,k))/(pplay(i,k))
4978            ELSE
4979               zx_qs = qsatl(t(i,k))/(pplay(i,k))
4980            ENDIF
4981            rh = q(i,k)/zx_qs
4982!
4983            d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
4984            d_rh = 1./tau*(rh_targ(i,k)-rh)
4985!
4986            tnew = t(i,k)+d_t(i,k)*dtime
4987!jyg<
4988!   Formule pour q :
4989!                         d_q = (1/tau) [rh_targ*qsat(T_new) - q] 
4990!
4991!  Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
4992!   qui n'etait pas correcte.
4993!
4994            IF (tnew.LT.RTT) THEN
4995               zx_qs = qsats(tnew)/(pplay(i,k))
4996            ELSE
4997               zx_qs = qsatl(tnew)/(pplay(i,k))
4998            ENDIF
4999!!            d_q(i,k) = d_q(i,k) + d_rh*zx_qs
5000            d_q(i,k) = d_q(i,k) + (1./tau)*(rh_targ(i,k)*zx_qs - q(i,k))
5001            rhnew = (q(i,k)+d_q(i,k)*dtime)/zx_qs
5002!
5003            print *,' k,d_t,rh,d_rh,rhnew,d_q ',    &
5004                      k,d_t(i,k),rh,d_rh,rhnew,d_q(i,k)
5005           ENDIF
5006!
5007         ENDDO
5008        ENDDO
5009!
5010      RETURN
5011      END
5012
5013      Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v,          &
5014     &                      d_u,d_v)
5015! ========================================================
5016  USE dimphy
5017
5018  implicit none
5019
5020! ========================================================
5021      REAL dtime
5022      REAL paprs(klon,klevp1)
5023      REAL pplay(klon,klev)
5024!
5025!      Variables d'etat
5026      REAL u(klon,klev)
5027      REAL v(klon,klev)
5028!
5029! Tendances
5030      REAL d_u(klon,klev)
5031      REAL d_v(klon,klev)
5032!
5033!   Profiles cible
5034      REAL u_targ(klon,klev)
5035      REAL v_targ(klon,klev)
5036!
5037!   Temps de relaxation
5038      REAL tau
5039!c      DATA tau /3600./
5040!      DATA tau /5400./
5041       DATA tau /43200./
5042!
5043   INTEGER k,i
5044
5045!
5046        print *,'dtime, tau ',dtime,tau
5047        print *, 'u_targ',u_targ
5048        print *, 'v_targ',v_targ
5049        print *,'zonal velocity ',u
5050        print *,'meridional velocity ',v
5051        DO k = 1,klev
5052         DO i = 1,klon
5053!CR: nudging everywhere
5054!           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
5055!
5056            d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
5057            d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
5058!
5059            print *,' k,u,d_u,v,d_v ',    &
5060                      k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
5061!           ENDIF
5062!
5063         ENDDO
5064        ENDDO
5065!
5066      RETURN
5067      END
5068
5069!=====================================================================
5070       SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas                                    &
5071     &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
5072     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
5073     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                   &
5074     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
5075     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
5076     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
5077!
5078     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
5079     &         ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
5080     &         ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                                          &
5081     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
5082     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
5083     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
5084 
5085       implicit none
5086 
5087#include "YOMCST.h"
5088#include "dimensions.h"
5089
5090!-------------------------------------------------------------------------
5091! Vertical interpolation of generic case forcing data onto mod_casel levels
5092!-------------------------------------------------------------------------
5093 
5094       integer nlevmax
5095       parameter (nlevmax=41)
5096       integer nlev_cas,mxcalc
5097!       real play(llm), plev_prof(nlevmax) 
5098!       real t_prof(nlevmax),q_prof(nlevmax)
5099!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
5100!       real ht_prof(nlevmax),vt_prof(nlevmax)
5101!       real hq_prof(nlevmax),vq_prof(nlevmax)
5102 
5103       real play(llm), plev_prof_cas(nlev_cas) 
5104       real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
5105       real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
5106       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
5107       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
5108       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
5109       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
5110       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
5111       real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
5112       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
5113 
5114       real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
5115       real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
5116       real u_mod_cas(llm),v_mod_cas(llm)
5117       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
5118       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
5119       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
5120       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
5121       real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
5122       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
5123 
5124       integer l,k,k1,k2
5125       real frac,frac1,frac2,fact
5126 
5127       do l = 1, llm
5128       print *,'debut interp2, play=',l,play(l)
5129       enddo
5130!      do l = 1, nlev_cas
5131!      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
5132!      enddo
5133
5134       do l = 1, llm
5135
5136        if (play(l).ge.plev_prof_cas(nlev_cas)) then
5137 
5138        mxcalc=l
5139        print *,'debut interp2, mxcalc=',mxcalc
5140         k1=0
5141         k2=0
5142
5143         if (play(l).le.plev_prof_cas(1)) then
5144
5145         do k = 1, nlev_cas-1
5146          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
5147            k1=k
5148            k2=k+1
5149          endif
5150         enddo
5151
5152         if (k1.eq.0 .or. k2.eq.0) then
5153          write(*,*) 'PB! k1, k2 = ',k1,k2
5154          write(*,*) 'l,play(l) = ',l,play(l)/100
5155         do k = 1, nlev_cas-1
5156          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
5157         enddo
5158         endif
5159
5160         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
5161         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
5162         theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
5163         t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
5164         thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
5165         thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
5166         qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
5167         ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
5168         qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
5169         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
5170         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
5171         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
5172         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
5173         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
5174         omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
5175         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
5176         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
5177         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
5178         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
5179         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
5180         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
5181         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
5182         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
5183         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
5184         dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
5185         hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
5186         vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
5187         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
5188         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
5189         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
5190         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
5191     
5192         else !play>plev_prof_cas(1)
5193
5194         k1=1
5195         k2=2
5196         print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
5197         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
5198         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
5199         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
5200         theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
5201         t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
5202         thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
5203         thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
5204         qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
5205         ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
5206         qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
5207         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
5208         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
5209         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
5210         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
5211         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
5212         omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
5213         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
5214         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
5215         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
5216         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
5217         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
5218         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
5219         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
5220         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
5221         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
5222         dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
5223         hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
5224         vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
5225         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
5226         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
5227         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
5228         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
5229
5230         endif ! play.le.plev_prof_cas(1)
5231
5232        else ! above max altitude of forcing file
5233 
5234!jyg
5235         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
5236         fact = max(fact,0.)                                           !jyg
5237         fact = exp(-fact)                                             !jyg
5238         t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
5239         theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
5240         thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
5241         thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
5242         qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
5243         ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
5244         qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
5245         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
5246         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
5247         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                     !jyg
5248         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                     !jyg
5249         w_mod_cas(l)= 0.0                                             !jyg
5250         omega_mod_cas(l)= 0.0                                         !jyg
5251         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
5252         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
5253         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
5254         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
5255         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
5256         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
5257         dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
5258         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
5259         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
5260         dth_mod_cas(l)= dth_prof_cas(nlev_cas) 
5261         hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
5262         vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
5263         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
5264         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
5265         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
5266         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
5267 
5268        endif ! play
5269 
5270       enddo ! l
5271
5272          return
5273          end
5274!***************************************************************************** 
5275
5276
Note: See TracBrowser for help on using the repository browser.