source: LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 3594

Last change on this file since 3594 was 3594, checked in by fhourdin, 5 years ago

Improvement for SCM standard format

  • Property svn:keywords set to Id
File size: 60.6 KB
Line 
1#include "conf_gcm.F90"
2
3!
4! $Id: 1DUTILS.h 3594 2019-10-29 06:32:03Z fhourdin $
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! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT
543!Config  Key  = nudging_q
544!Config  Desc = forcage ou non par nudging sur q
545!Config  Def  = false
546!Config  Help = forcage ou non par nudging sur q
547       nudging_qv =0
548       CALL getin('nudging_q',nudging_qv)
549       CALL getin('nudging_qv',nudging_qv)
550
551       p_nudging_u=11000.
552       p_nudging_v=11000.
553       p_nudging_t=11000.
554       p_nudging_qv=11000.
555       CALL getin('p_nudging_u',p_nudging_u)
556       CALL getin('p_nudging_v',p_nudging_v)
557       CALL getin('p_nudging_t',p_nudging_t)
558       CALL getin('p_nudging_qv',p_nudging_qv)
559
560!Config  Key  = nudging_t
561!Config  Desc = forcage ou non par nudging sur t
562!Config  Def  = false
563!Config  Help = forcage ou non par nudging sur t
564       nudging_t =0
565       CALL getin('nudging_t',nudging_t)
566
567
568
569      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
570      write(lunout,*)' Configuration des parametres du gcm1D: '
571      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
572      write(lunout,*)' restart = ', restart
573      write(lunout,*)' forcing_type = ', forcing_type
574      write(lunout,*)' time_ini = ', time_ini
575      write(lunout,*)' rlat = ', xlat
576      write(lunout,*)' rlon = ', xlon
577      write(lunout,*)' airephy = ', airefi
578      write(lunout,*)' nat_surf = ', nat_surf
579      write(lunout,*)' tsurf = ', tsurf
580      write(lunout,*)' psurf = ', psurf
581      write(lunout,*)' zsurf = ', zsurf
582      write(lunout,*)' rugos = ', rugos
583      write(lunout,*)' snowmass=', snowmass
584      write(lunout,*)' wtsurf = ', wtsurf
585      write(lunout,*)' wqsurf = ', wqsurf
586      write(lunout,*)' albedo = ', albedo
587      write(lunout,*)' xagesno = ', xagesno
588      write(lunout,*)' restart_runoff = ', restart_runoff
589      write(lunout,*)' qsolinp = ', qsolinp
590      write(lunout,*)' zpicinp = ', zpicinp
591      write(lunout,*)' nudge_tsoil = ', nudge_tsoil
592      write(lunout,*)' isoil_nudge = ', isoil_nudge
593      write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
594      write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
595      write(lunout,*)' tadv =      ', tadv
596      write(lunout,*)' tadvv =     ', tadvv
597      write(lunout,*)' tadvh =     ', tadvh
598      write(lunout,*)' thadv =     ', thadv
599      write(lunout,*)' thadvv =    ', thadvv
600      write(lunout,*)' thadvh =    ', thadvh
601      write(lunout,*)' qadv =      ', qadv
602      write(lunout,*)' qadvv =     ', qadvv
603      write(lunout,*)' qadvh =     ', qadvh
604      write(lunout,*)' trad =      ', trad
605      write(lunout,*)' forc_omega = ', forc_omega
606      write(lunout,*)' forc_w     = ', forc_w
607      write(lunout,*)' forc_geo   = ', forc_geo
608      write(lunout,*)' forc_ustar = ', forc_ustar
609      write(lunout,*)' nudging_u  = ', nudging_u
610      write(lunout,*)' nudging_v  = ', nudging_v
611      write(lunout,*)' nudging_t  = ', nudging_t
612      write(lunout,*)' nudging_qv  = ', nudging_qv
613      IF (forcing_type .eq.40) THEN
614        write(lunout,*) '--- Forcing type GCSS Old --- with:'
615        write(lunout,*)'imp_fcg',imp_fcg_gcssold
616        write(lunout,*)'ts_fcg',ts_fcg_gcssold
617        write(lunout,*)'tp_fcg',Tp_fcg_gcssold
618        write(lunout,*)'tp_ini',Tp_ini_gcssold
619        write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
620      ENDIF
621
622      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
623      write(lunout,*)
624!
625      RETURN
626      END
627!
628! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
629!
630!
631      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,                 &
632     &                          ucov,vcov,temp,q,omega2)
633      USE dimphy
634      USE mod_grid_phy_lmdz
635      USE mod_phys_lmdz_para
636      USE iophy
637      USE phys_state_var_mod
638      USE iostart
639      USE write_field_phy
640      USE infotrac
641      use control_mod
642      USE comconst_mod, ONLY: im, jm, lllm
643      USE logic_mod, ONLY: fxyhypb, ysinus
644      USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn
645
646      IMPLICIT NONE
647!=======================================================
648! Ecriture du fichier de redemarrage sous format NetCDF
649!=======================================================
650!   Declarations:
651!   -------------
652      include "dimensions.h"
653!!#include "control.h"
654      include "netcdf.inc"
655
656!   Arguments:
657!   ----------
658      CHARACTER*(*) fichnom
659!Al1 plev tronque pour .nc mais plev(klev+1):=0
660      real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev)
661      real :: presnivs(klon,klev)
662      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
663      real :: q(klon,klev,nqtot),omega2(klon,klev)
664!      real :: ug(klev),vg(klev),fcoriolis
665      real :: phis(klon) 
666
667!   Variables locales pour NetCDF:
668!   ------------------------------
669      INTEGER iq
670      INTEGER length
671      PARAMETER (length = 100)
672      REAL tab_cntrl(length) ! tableau des parametres du run
673      character*4 nmq(nqtot)
674      character*12 modname
675      character*80 abort_message
676      LOGICAL found
677
678      modname = 'dyn1deta0 : '
679!!      nmq(1)="vap"
680!!      nmq(2)="cond"
681!!      do iq=3,nqtot
682!!        write(nmq(iq),'("tra",i1)') iq-2
683!!      enddo
684      DO iq = 1,nqtot
685        nmq(iq) = trim(tname(iq))
686      ENDDO
687      print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
688      CALL open_startphy(fichnom)
689      print*,'after open startphy ',fichnom,nmq
690
691!
692! Lecture des parametres de controle:
693!
694      CALL get_var("controle",tab_cntrl)
695       
696
697      im         = tab_cntrl(1)
698      jm         = tab_cntrl(2)
699      lllm       = tab_cntrl(3)
700      day_ref    = tab_cntrl(4)
701      annee_ref  = tab_cntrl(5)
702!      rad        = tab_cntrl(6)
703!      omeg       = tab_cntrl(7)
704!      g          = tab_cntrl(8)
705!      cpp        = tab_cntrl(9)
706!      kappa      = tab_cntrl(10)
707!      daysec     = tab_cntrl(11)
708!      dtvr       = tab_cntrl(12)
709!      etot0      = tab_cntrl(13)
710!      ptot0      = tab_cntrl(14)
711!      ztot0      = tab_cntrl(15)
712!      stot0      = tab_cntrl(16)
713!      ang0       = tab_cntrl(17)
714!      pa         = tab_cntrl(18)
715!      preff      = tab_cntrl(19)
716!
717!      clon       = tab_cntrl(20)
718!      clat       = tab_cntrl(21)
719!      grossismx  = tab_cntrl(22)
720!      grossismy  = tab_cntrl(23)
721!
722      IF ( tab_cntrl(24).EQ.1. )  THEN
723        fxyhypb  =.true.
724!        dzoomx   = tab_cntrl(25)
725!        dzoomy   = tab_cntrl(26)
726!        taux     = tab_cntrl(28)
727!        tauy     = tab_cntrl(29)
728      ELSE
729        fxyhypb = .false.
730        ysinus  = .false.
731        IF( tab_cntrl(27).EQ.1. ) ysinus =.true. 
732      ENDIF
733
734      day_ini = tab_cntrl(30)
735      itau_dyn = tab_cntrl(31)
736!   .................................................................
737!
738!
739!      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
740!Al1
741       Print*,'day_ref,annee_ref,day_ini,itau_dyn',                         &
742     &              day_ref,annee_ref,day_ini,itau_dyn
743
744!  Lecture des champs
745!
746      CALL get_field("play",play,found)
747      IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
748      CALL get_field("phi",phi,found)
749      IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
750      CALL get_field("phis",phis,found)
751      IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
752      CALL get_field("presnivs",presnivs,found)
753      IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
754      CALL get_field("ucov",ucov,found)
755      IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
756      CALL get_field("vcov",vcov,found)
757      IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
758      CALL get_field("temp",temp,found)
759      IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
760      CALL get_field("omega2",omega2,found)
761      IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
762      plev(1,klev+1)=0.
763      CALL get_field("plev",plev(:,1:klev),found)
764      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
765
766      Do iq=1,nqtot
767        CALL get_field("q"//nmq(iq),q(:,:,iq),found)
768        IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
769      EndDo
770
771      CALL close_startphy
772      print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
773!
774      RETURN
775      END
776!
777! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
778!
779!
780      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,           &
781     &                          ucov,vcov,temp,q,omega2)
782      USE dimphy
783      USE mod_grid_phy_lmdz
784      USE mod_phys_lmdz_para
785      USE phys_state_var_mod
786      USE iostart
787      USE infotrac
788      use control_mod
789      USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
790      USE logic_mod, ONLY: fxyhypb, ysinus
791      USE temps_mod, ONLY: annee_ref,day_end,day_ref,itau_dyn,itaufin
792
793      IMPLICIT NONE
794!=======================================================
795! Ecriture du fichier de redemarrage sous format NetCDF
796!=======================================================
797!   Declarations:
798!   -------------
799      include "dimensions.h"
800!!#include "control.h"
801      include "netcdf.inc"
802
803!   Arguments:
804!   ----------
805      CHARACTER*(*) fichnom
806!Al1 plev tronque pour .nc mais plev(klev+1):=0
807      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
808      real :: presnivs(klon,klev)
809      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
810      real :: q(klon,klev,nqtot)
811      real :: omega2(klon,klev),rho(klon,klev+1)
812!      real :: ug(klev),vg(klev),fcoriolis
813      real :: phis(klon) 
814
815!   Variables locales pour NetCDF:
816!   ------------------------------
817      INTEGER nid
818      INTEGER ierr
819      INTEGER iq,l
820      INTEGER length
821      PARAMETER (length = 100)
822      REAL tab_cntrl(length) ! tableau des parametres du run
823      character*4 nmq(nqtot)
824      character*20 modname
825      character*80 abort_message
826!
827      INTEGER pass
828
829      CALL open_restartphy(fichnom)
830      print*,'redm1 ',fichnom,klon,klev,nqtot
831!!      nmq(1)="vap"
832!!      nmq(2)="cond"
833!!      nmq(3)="tra1"
834!!      nmq(4)="tra2"
835      DO iq = 1,nqtot
836        nmq(iq) = trim(tname(iq))
837      ENDDO
838
839!     modname = 'dyn1dredem'
840!     ierr = NF_OPEN(fichnom, NF_WRITE, nid)
841!     IF (ierr .NE. NF_NOERR) THEN
842!        abort_message="Pb. d ouverture "//fichnom
843!        CALL abort_gcm('Modele 1D',abort_message,1)
844!     ENDIF
845
846      DO l=1,length
847       tab_cntrl(l) = 0.
848      ENDDO
849       tab_cntrl(1)  = FLOAT(iim)
850       tab_cntrl(2)  = FLOAT(jjm)
851       tab_cntrl(3)  = FLOAT(llm)
852       tab_cntrl(4)  = FLOAT(day_ref)
853       tab_cntrl(5)  = FLOAT(annee_ref)
854       tab_cntrl(6)  = rad
855       tab_cntrl(7)  = omeg
856       tab_cntrl(8)  = g
857       tab_cntrl(9)  = cpp
858       tab_cntrl(10) = kappa
859       tab_cntrl(11) = daysec
860       tab_cntrl(12) = dtvr
861!       tab_cntrl(13) = etot0
862!       tab_cntrl(14) = ptot0
863!       tab_cntrl(15) = ztot0
864!       tab_cntrl(16) = stot0
865!       tab_cntrl(17) = ang0
866!       tab_cntrl(18) = pa
867!       tab_cntrl(19) = preff
868!
869!    .....    parametres  pour le zoom      ......   
870
871!       tab_cntrl(20)  = clon
872!       tab_cntrl(21)  = clat
873!       tab_cntrl(22)  = grossismx
874!       tab_cntrl(23)  = grossismy
875!
876      IF ( fxyhypb )   THEN
877       tab_cntrl(24) = 1.
878!       tab_cntrl(25) = dzoomx
879!       tab_cntrl(26) = dzoomy
880       tab_cntrl(27) = 0.
881!       tab_cntrl(28) = taux
882!       tab_cntrl(29) = tauy
883      ELSE
884       tab_cntrl(24) = 0.
885!       tab_cntrl(25) = dzoomx
886!       tab_cntrl(26) = dzoomy
887       tab_cntrl(27) = 0.
888       tab_cntrl(28) = 0.
889       tab_cntrl(29) = 0.
890       IF( ysinus )  tab_cntrl(27) = 1.
891      ENDIF
892!Al1 iday_end -> day_end
893       tab_cntrl(30) = FLOAT(day_end)
894       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
895!
896      DO pass=1,2
897      CALL put_var(pass,"controle","Param. de controle Dyn1D",tab_cntrl)
898!
899
900!  Ecriture/extension de la coordonnee temps
901
902
903!  Ecriture des champs
904!
905      CALL put_field(pass,"plev","p interfaces sauf la nulle",plev)
906      CALL put_field(pass,"play","",play)
907      CALL put_field(pass,"phi","geopotentielle",phi)
908      CALL put_field(pass,"phis","geopotentiell de surface",phis)
909      CALL put_field(pass,"presnivs","",presnivs)
910      CALL put_field(pass,"ucov","",ucov)
911      CALL put_field(pass,"vcov","",vcov)
912      CALL put_field(pass,"temp","",temp)
913      CALL put_field(pass,"omega2","",omega2)
914
915      Do iq=1,nqtot
916        CALL put_field(pass,"q"//nmq(iq),"eau vap ou condens et traceurs",           &
917     &                                                      q(:,:,iq))
918      EndDo
919    IF (pass==1) CALL enddef_restartphy
920    IF (pass==2) CALL close_restartphy
921
922
923      ENDDO
924
925!
926      RETURN
927      END
928      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
929      IMPLICIT NONE
930!=======================================================================
931!   passage d'un champ de la grille scalaire a la grille physique
932!=======================================================================
933 
934!-----------------------------------------------------------------------
935!   declarations:
936!   -------------
937 
938      INTEGER im,jm,ngrid,nfield
939      REAL pdyn(im,jm,nfield)
940      REAL pfi(ngrid,nfield)
941 
942      INTEGER i,j,ifield,ig
943 
944!-----------------------------------------------------------------------
945!   calcul:
946!   -------
947 
948      DO ifield=1,nfield
949!   traitement des poles
950         DO i=1,im
951            pdyn(i,1,ifield)=pfi(1,ifield)
952            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
953         ENDDO
954 
955!   traitement des point normaux
956         DO j=2,jm-1
957            ig=2+(j-2)*(im-1)
958            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
959            pdyn(im,j,ifield)=pdyn(1,j,ifield)
960         ENDDO
961      ENDDO
962 
963      RETURN
964      END
965 
966 
967
968      SUBROUTINE abort_gcm(modname, message, ierr)
969 
970      USE IOIPSL
971!
972! Stops the simulation cleanly, closing files and printing various
973! comments
974!
975!  Input: modname = name of calling program
976!         message = stuff to print
977!         ierr    = severity of situation ( = 0 normal )
978 
979      character(len=*) modname
980      integer ierr
981      character(len=*) message
982 
983      write(*,*) 'in abort_gcm'
984      call histclo
985!     call histclo(2)
986!     call histclo(3)
987!     call histclo(4)
988!     call histclo(5)
989      write(*,*) 'out of histclo'
990      write(*,*) 'Stopping in ', modname
991      write(*,*) 'Reason = ',message
992      call getin_dump
993!
994      if (ierr .eq. 0) then
995        write(*,*) 'Everything is cool'
996      else
997        write(*,*) 'Houston, we have a problem ', ierr
998      endif
999      STOP
1000      END
1001      REAL FUNCTION fq_sat(kelvin, millibar)
1002!
1003      IMPLICIT none
1004!======================================================================
1005! Autheur(s): Z.X. Li (LMD/CNRS)
1006! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
1007!======================================================================
1008! Arguments:
1009! kelvin---input-R: temperature en Kelvin
1010! millibar--input-R: pression en mb
1011!
1012! fq_sat----output-R: vapeur d'eau saturante en kg/kg
1013!======================================================================
1014!
1015      REAL kelvin, millibar
1016!
1017      REAL r2es
1018      PARAMETER (r2es=611.14 *18.0153/28.9644)
1019!
1020      REAL r3les, r3ies, r3es
1021      PARAMETER (R3LES=17.269)
1022      PARAMETER (R3IES=21.875)
1023!
1024      REAL r4les, r4ies, r4es
1025      PARAMETER (R4LES=35.86)
1026      PARAMETER (R4IES=7.66)
1027!
1028      REAL rtt
1029      PARAMETER (rtt=273.16)
1030!
1031      REAL retv
1032      PARAMETER (retv=28.9644/18.0153 - 1.0)
1033!
1034      REAL zqsat
1035      REAL temp, pres
1036!     ------------------------------------------------------------------
1037!
1038!
1039      temp = kelvin
1040      pres = millibar * 100.0
1041!      write(*,*)'kelvin,millibar=',kelvin,millibar
1042!      write(*,*)'temp,pres=',temp,pres
1043!
1044      IF (temp .LE. rtt) THEN
1045         r3es = r3ies
1046         r4es = r4ies
1047      ELSE
1048         r3es = r3les
1049         r4es = r4les
1050      ENDIF
1051!
1052      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
1053      zqsat=MIN(0.5,ZQSAT)
1054      zqsat=zqsat/(1.-retv  *zqsat)
1055!
1056      fq_sat = zqsat
1057!
1058      RETURN
1059      END
1060 
1061      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
1062      IMPLICIT NONE
1063!=======================================================================
1064!   passage d'un champ de la grille scalaire a la grille physique
1065!=======================================================================
1066 
1067!-----------------------------------------------------------------------
1068!   declarations:
1069!   -------------
1070 
1071      INTEGER im,jm,ngrid,nfield
1072      REAL pdyn(im,jm,nfield)
1073      REAL pfi(ngrid,nfield)
1074 
1075      INTEGER j,ifield,ig
1076 
1077!-----------------------------------------------------------------------
1078!   calcul:
1079!   -------
1080 
1081      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)                          &
1082     &    STOP 'probleme de dim'
1083!   traitement des poles
1084      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
1085      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
1086 
1087!   traitement des point normaux
1088      DO ifield=1,nfield
1089         DO j=2,jm-1
1090            ig=2+(j-2)*(im-1)
1091            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
1092         ENDDO
1093      ENDDO
1094 
1095      RETURN
1096      END
1097 
1098      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
1099 
1100!    Ancienne version disvert dont on a modifie nom pour utiliser
1101!    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
1102!    (MPL 18092012)
1103!
1104!    Auteur :  P. Le Van .
1105!
1106      IMPLICIT NONE
1107 
1108      include "dimensions.h"
1109      include "paramet.h"
1110!
1111!=======================================================================
1112!
1113!
1114!    s = sigma ** kappa   :  coordonnee  verticale
1115!    dsig(l)            : epaisseur de la couche l ds la coord.  s
1116!    sig(l)             : sigma a l'interface des couches l et l-1
1117!    ds(l)              : distance entre les couches l et l-1 en coord.s
1118!
1119!=======================================================================
1120!
1121      REAL pa,preff
1122      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
1123      REAL presnivs(llm)
1124!
1125!   declarations:
1126!   -------------
1127!
1128      REAL sig(llm+1),dsig(llm)
1129!
1130      INTEGER l
1131      REAL snorm
1132      REAL alpha,beta,gama,delta,deltaz,h
1133      INTEGER np,ierr
1134      REAL pi,x
1135 
1136!-----------------------------------------------------------------------
1137!
1138      pi=2.*ASIN(1.)
1139 
1140      OPEN(99,file='sigma.def',status='old',form='formatted',                   &
1141     &   iostat=ierr)
1142 
1143!-----------------------------------------------------------------------
1144!   cas 1 on lit les options dans sigma.def:
1145!   ----------------------------------------
1146 
1147      IF (ierr.eq.0) THEN
1148 
1149      print*,'WARNING!!! on lit les options dans sigma.def'
1150      READ(99,*) deltaz
1151      READ(99,*) h
1152      READ(99,*) beta
1153      READ(99,*) gama
1154      READ(99,*) delta
1155      READ(99,*) np
1156      CLOSE(99)
1157      alpha=deltaz/(llm*h)
1158!
1159 
1160       DO 1  l = 1, llm
1161       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*                    &
1162     &          ( (tanh(gama*l)/tanh(gama*llm))**np +                      &
1163     &            (1.-l/FLOAT(llm))*delta )
1164   1   CONTINUE
1165 
1166       sig(1)=1.
1167       DO 101 l=1,llm-1
1168          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
1169101    CONTINUE
1170       sig(llm+1)=0.
1171 
1172       DO 2  l = 1, llm
1173       dsig(l) = sig(l)-sig(l+1)
1174   2   CONTINUE
1175!
1176 
1177      ELSE
1178!-----------------------------------------------------------------------
1179!   cas 2 ancienne discretisation (LMD5...):
1180!   ----------------------------------------
1181 
1182      PRINT*,'WARNING!!! Ancienne discretisation verticale'
1183 
1184      h=7.
1185      snorm  = 0.
1186      DO l = 1, llm
1187         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
1188         dsig(l) = 1.0 + 7.0 * SIN(x)**2
1189         snorm = snorm + dsig(l)
1190      ENDDO
1191      snorm = 1./snorm
1192      DO l = 1, llm
1193         dsig(l) = dsig(l)*snorm
1194      ENDDO
1195      sig(llm+1) = 0.
1196      DO l = llm, 1, -1
1197         sig(l) = sig(l+1) + dsig(l)
1198      ENDDO
1199 
1200      ENDIF
1201 
1202 
1203      DO l=1,llm
1204        nivsigs(l) = FLOAT(l)
1205      ENDDO
1206 
1207      DO l=1,llmp1
1208        nivsig(l)= FLOAT(l)
1209      ENDDO
1210 
1211!
1212!    ....  Calculs  de ap(l) et de bp(l)  ....
1213!    .........................................
1214!
1215!
1216!   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
1217!
1218 
1219      bp(llmp1) =   0.
1220 
1221      DO l = 1, llm
1222!c
1223!cc    ap(l) = 0.
1224!cc    bp(l) = sig(l)
1225 
1226      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
1227      ap(l) = pa * ( sig(l) - bp(l) )
1228!
1229      ENDDO
1230      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
1231 
1232      PRINT *,' BP '
1233      PRINT *,  bp
1234      PRINT *,' AP '
1235      PRINT *,  ap
1236 
1237      DO l = 1, llm
1238       dpres(l) = bp(l) - bp(l+1)
1239       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
1240      ENDDO
1241 
1242      PRINT *,' PRESNIVS '
1243      PRINT *,presnivs
1244 
1245      RETURN
1246      END
1247
1248!======================================================================
1249       SUBROUTINE read_tsurf1d(knon,sst_out)
1250
1251! This subroutine specifies the surface temperature to be used in 1D simulations
1252
1253      USE dimphy, ONLY : klon
1254
1255      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
1256      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
1257
1258       INTEGER :: i
1259! COMMON defined in lmdz1d.F:
1260       real ts_cur
1261       common /sst_forcing/ts_cur
1262
1263       DO i = 1, knon
1264        sst_out(i) = ts_cur
1265       ENDDO
1266
1267      END SUBROUTINE read_tsurf1d
1268
1269!===============================================================
1270      subroutine advect_vert(llm,w,dt,q,plev)
1271!===============================================================
1272!   Schema amont pour l'advection verticale en 1D
1273!   w est la vitesse verticale dp/dt en Pa/s
1274!   Traitement en volumes finis
1275!   d / dt ( zm q ) = delta_z ( omega q )
1276!   d / dt ( zm ) = delta_z ( omega )
1277!   avec zm = delta_z ( p )
1278!   si * designe la valeur au pas de temps t+dt
1279!   zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1280!   zm*(l) -zm(l) = w(l+1) - w(l)
1281!   avec w=omega * dt
1282!---------------------------------------------------------------
1283      implicit none
1284! arguments
1285      integer llm
1286      real w(llm+1),q(llm),plev(llm+1),dt
1287
1288! local
1289      integer l
1290      real zwq(llm+1),zm(llm+1),zw(llm+1)
1291      real qold
1292
1293!---------------------------------------------------------------
1294
1295      do l=1,llm
1296         zw(l)=dt*w(l)
1297         zm(l)=plev(l)-plev(l+1)
1298         zwq(l)=q(l)*zw(l)
1299      enddo
1300      zwq(llm+1)=0.
1301      zw(llm+1)=0.
1302 
1303      do l=1,llm
1304         qold=q(l)
1305         q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1306         print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1307      enddo
1308
1309 
1310      return
1311      end
1312
1313!===============================================================
1314
1315
1316       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,              &
1317     &                q,temp,u,v,play)
1318!itlmd
1319!----------------------------------------------------------------------
1320!   Calcul de l'advection verticale (ascendance et subsidence) de
1321!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1322!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1323!   sans WTG rajouter une advection horizontale
1324!---------------------------------------------------------------------- 
1325        implicit none
1326#include "YOMCST.h"
1327!        argument
1328        integer llm
1329        real  omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1330        real  d_u_va(llm), d_v_va(llm)
1331        real  q(llm,3),temp(llm)
1332        real  u(llm),v(llm)
1333        real  play(llm)
1334! interne
1335        integer l
1336        real alpha,omgdown,omgup
1337
1338      do l= 1,llm
1339       if(l.eq.1) then
1340!si omgup pour la couche 1, alors tendance nulle
1341        omgdown=max(omega(2),0.0)
1342        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1343        d_t_va(l)= alpha*(omgdown)-omgdown*(temp(l)-temp(l+1))             &
1344     &       /(play(l)-play(l+1))
1345
1346        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))             
1347
1348        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))             
1349        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))             
1350
1351       
1352       elseif(l.eq.llm) then
1353        omgup=min(omega(l),0.0)
1354        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1355        d_t_va(l)= alpha*(omgup)-                                          &
1356
1357!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1358     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1359        d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1360        d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1361        d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1362       
1363       else
1364        omgup=min(omega(l),0.0)
1365        omgdown=max(omega(l+1),0.0)
1366        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1367        d_t_va(l)= alpha*(omgup+omgdown)-omgdown*(temp(l)-temp(l+1))       &
1368     &              /(play(l)-play(l+1))-                                  &
1369!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1370     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1371!      print*, '  ??? '
1372
1373        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))                            &
1374     &              /(play(l)-play(l+1))-                                  &
1375     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) 
1376        d_u_va(l)= -omgdown*(u(l)-u(l+1))                                  &
1377     &              /(play(l)-play(l+1))-                                  &
1378     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) 
1379        d_v_va(l)= -omgdown*(v(l)-v(l+1))                                  &
1380     &              /(play(l)-play(l+1))-                                  &
1381     &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1382       
1383      endif
1384         
1385      enddo
1386!fin itlmd
1387        return
1388        end
1389!       SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1390       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,                        &
1391     &                q,temp,u,v,play)
1392!itlmd
1393!----------------------------------------------------------------------
1394!   Calcul de l'advection verticale (ascendance et subsidence) de
1395!   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1396!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1397!   sans WTG rajouter une advection horizontale
1398!---------------------------------------------------------------------- 
1399        implicit none
1400#include "YOMCST.h"
1401!        argument
1402        integer llm,nqtot
1403        real  omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1404!        real  d_u_va(llm), d_v_va(llm)
1405        real  q(llm,nqtot),temp(llm)
1406        real  u(llm),v(llm)
1407        real  play(llm)
1408        real cor(llm)
1409!        real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1410        real dph(llm),dqdp(llm),dtdp(llm)
1411! interne
1412        integer k
1413        real omdn,omup
1414
1415!        dudp=0.
1416!        dvdp=0.
1417        dqdp=0.
1418        dtdp=0.
1419!        d_u_va=0.
1420!        d_v_va=0.
1421
1422      cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
1423
1424
1425      do k=2,llm-1
1426
1427       dph  (k-1) = (play()- play(k-1  ))
1428!       dudp (k-1) = (u   ()- u   (k-1  ))/dph(k-1)
1429!       dvdp (k-1) = (v   ()- v   (k-1  ))/dph(k-1)
1430       dqdp (k-1) = (q   (k,1)- q   (k-1,1))/dph(k-1)
1431       dtdp (k-1) = (temp()- temp(k-1  ))/dph(k-1)
1432
1433      enddo
1434
1435!      dudp (  llm  ) = dudp ( llm-1 )
1436!      dvdp (  llm  ) = dvdp ( llm-1 )
1437      dqdp (  llm  ) = dqdp ( llm-1 )
1438      dtdp (  llm  ) = dtdp ( llm-1 )
1439
1440      do k=2,llm-1
1441      omdn=max(0.0,omega(k+1))
1442      omup=min(0.0,omega( k ))
1443
1444!      d_u_va(k)  = -omdn*dudp(k)-omup*dudp(k-1)
1445!      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
1446      d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1)
1447      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)+(omup+omdn)*cor(k)
1448      enddo
1449
1450      omdn=max(0.0,omega( 2 ))
1451      omup=min(0.0,omega(llm))
1452!      d_u_va( 1 )   = -omdn*dudp( 1 )
1453!      d_u_va(llm)   = -omup*dudp(llm)
1454!      d_v_va( 1 )   = -omdn*dvdp( 1 )
1455!      d_v_va(llm)   = -omup*dvdp(llm)
1456      d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1457      d_q_va(llm,1) = -omup*dqdp(llm)
1458      d_t_va( 1 )   = -omdn*dtdp( 1 )+omdn*cor( 1 )
1459      d_t_va(llm)   = -omup*dtdp(llm)!+omup*cor(llm)
1460
1461!      if(abs(rlat(1))>10.) then
1462!     Calculate the tendency due agestrophic motions
1463!      du_age = fcoriolis*(v-vg)
1464!      dv_age = fcoriolis*(ug-u)
1465!      endif
1466
1467!       call writefield_phy('d_t_va',d_t_va,llm)
1468
1469          return
1470         end
1471
1472!======================================================================
1473
1474!  Subroutines for nudging
1475
1476      Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
1477! ========================================================
1478  USE dimphy
1479
1480  implicit none
1481
1482! ========================================================
1483      REAL paprs(klon,klevp1)
1484      REAL pplay(klon,klev)
1485!
1486!      Variables d'etat
1487      REAL t(klon,klev)
1488      REAL q(klon,klev)
1489!
1490!   Profiles cible
1491      REAL t_targ(klon,klev)
1492      REAL rh_targ(klon,klev)
1493!
1494   INTEGER k,i
1495   REAL zx_qs
1496
1497! Declaration des constantes et des fonctions thermodynamiques
1498!
1499include "YOMCST.h"
1500include "YOETHF.h"
1501!
1502!  ----------------------------------------
1503!  Statement functions
1504include "FCTTRE.h"
1505!  ----------------------------------------
1506!
1507        DO k = 1,klev
1508         DO i = 1,klon
1509           t_targ(i,k) = t(i,k)
1510           IF (t(i,k).LT.RTT) THEN
1511              zx_qs = qsats(t(i,k))/(pplay(i,k))
1512           ELSE
1513              zx_qs = qsatl(t(i,k))/(pplay(i,k))
1514           ENDIF
1515           rh_targ(i,k) = q(i,k)/zx_qs
1516         ENDDO
1517        ENDDO
1518      print *, 't_targ',t_targ
1519      print *, 'rh_targ',rh_targ
1520!
1521!
1522      RETURN
1523      END
1524
1525      Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
1526! ========================================================
1527  USE dimphy
1528
1529  implicit none
1530
1531! ========================================================
1532      REAL paprs(klon,klevp1)
1533      REAL pplay(klon,klev)
1534!
1535!      Variables d'etat
1536      REAL u(klon,klev)
1537      REAL v(klon,klev)
1538!
1539!   Profiles cible
1540      REAL u_targ(klon,klev)
1541      REAL v_targ(klon,klev)
1542!
1543   INTEGER k,i
1544!
1545        DO k = 1,klev
1546         DO i = 1,klon
1547           u_targ(i,k) = u(i,k)
1548           v_targ(i,k) = v(i,k)
1549         ENDDO
1550        ENDDO
1551      print *, 'u_targ',u_targ
1552      print *, 'v_targ',v_targ
1553!
1554!
1555      RETURN
1556      END
1557
1558      Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q,          &
1559     &                      d_t,d_q)
1560! ========================================================
1561  USE dimphy
1562
1563  implicit none
1564
1565! ========================================================
1566      REAL dtime
1567      REAL paprs(klon,klevp1)
1568      REAL pplay(klon,klev)
1569!
1570!      Variables d'etat
1571      REAL t(klon,klev)
1572      REAL q(klon,klev)
1573!
1574! Tendances
1575      REAL d_t(klon,klev)
1576      REAL d_q(klon,klev)
1577!
1578!   Profiles cible
1579      REAL t_targ(klon,klev)
1580      REAL rh_targ(klon,klev)
1581!
1582!   Temps de relaxation
1583      REAL tau
1584!c      DATA tau /3600./
1585!!      DATA tau /5400./
1586      DATA tau /1800./
1587!
1588   INTEGER k,i
1589   REAL zx_qs, rh, tnew, d_rh, rhnew
1590
1591! Declaration des constantes et des fonctions thermodynamiques
1592!
1593include "YOMCST.h"
1594include "YOETHF.h"
1595!
1596!  ----------------------------------------
1597!  Statement functions
1598include "FCTTRE.h"
1599!  ----------------------------------------
1600!
1601        print *,'dtime, tau ',dtime,tau
1602        print *, 't_targ',t_targ
1603        print *, 'rh_targ',rh_targ
1604        print *,'temp ',t
1605        print *,'hum ',q
1606!
1607        DO k = 1,klev
1608         DO i = 1,klon
1609           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
1610            IF (t(i,k).LT.RTT) THEN
1611               zx_qs = qsats(t(i,k))/(pplay(i,k))
1612            ELSE
1613               zx_qs = qsatl(t(i,k))/(pplay(i,k))
1614            ENDIF
1615            rh = q(i,k)/zx_qs
1616!
1617            d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
1618            d_rh = 1./tau*(rh_targ(i,k)-rh)
1619!
1620            tnew = t(i,k)+d_t(i,k)*dtime
1621!jyg<
1622!   Formule pour q :
1623!                         d_q = (1/tau) [rh_targ*qsat(T_new) - q] 
1624!
1625!  Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
1626!   qui n'etait pas correcte.
1627!
1628            IF (tnew.LT.RTT) THEN
1629               zx_qs = qsats(tnew)/(pplay(i,k))
1630            ELSE
1631               zx_qs = qsatl(tnew)/(pplay(i,k))
1632            ENDIF
1633!!            d_q(i,k) = d_q(i,k) + d_rh*zx_qs
1634            d_q(i,k) = d_q(i,k) + (1./tau)*(rh_targ(i,k)*zx_qs - q(i,k))
1635            rhnew = (q(i,k)+d_q(i,k)*dtime)/zx_qs
1636!
1637            print *,' k,d_t,rh,d_rh,rhnew,d_q ',    &
1638                      k,d_t(i,k),rh,d_rh,rhnew,d_q(i,k)
1639           ENDIF
1640!
1641         ENDDO
1642        ENDDO
1643!
1644      RETURN
1645      END
1646
1647      Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v,          &
1648     &                      d_u,d_v)
1649! ========================================================
1650  USE dimphy
1651
1652  implicit none
1653
1654! ========================================================
1655      REAL dtime
1656      REAL paprs(klon,klevp1)
1657      REAL pplay(klon,klev)
1658!
1659!      Variables d'etat
1660      REAL u(klon,klev)
1661      REAL v(klon,klev)
1662!
1663! Tendances
1664      REAL d_u(klon,klev)
1665      REAL d_v(klon,klev)
1666!
1667!   Profiles cible
1668      REAL u_targ(klon,klev)
1669      REAL v_targ(klon,klev)
1670!
1671!   Temps de relaxation
1672      REAL tau
1673!c      DATA tau /3600./
1674!      DATA tau /5400./
1675       DATA tau /43200./
1676!
1677   INTEGER k,i
1678
1679!
1680        print *,'dtime, tau ',dtime,tau
1681        print *, 'u_targ',u_targ
1682        print *, 'v_targ',v_targ
1683        print *,'zonal velocity ',u
1684        print *,'meridional velocity ',v
1685        DO k = 1,klev
1686         DO i = 1,klon
1687!CR: nudging everywhere
1688!           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
1689!
1690            d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
1691            d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
1692!
1693            print *,' k,u,d_u,v,d_v ',    &
1694                      k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
1695!           ENDIF
1696!
1697         ENDDO
1698        ENDDO
1699!
1700      RETURN
1701      END
1702
1703!=====================================================================
1704       SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas                                    &
1705     &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
1706     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
1707     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                   &
1708     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
1709     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
1710     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
1711!
1712     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
1713     &         ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
1714     &         ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                                          &
1715     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
1716     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
1717     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
1718 
1719       implicit none
1720 
1721#include "YOMCST.h"
1722#include "dimensions.h"
1723
1724!-------------------------------------------------------------------------
1725! Vertical interpolation of generic case forcing data onto mod_casel levels
1726!-------------------------------------------------------------------------
1727 
1728       integer nlevmax
1729       parameter (nlevmax=41)
1730       integer nlev_cas,mxcalc
1731!       real play(llm), plev_prof(nlevmax) 
1732!       real t_prof(nlevmax),q_prof(nlevmax)
1733!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1734!       real ht_prof(nlevmax),vt_prof(nlevmax)
1735!       real hq_prof(nlevmax),vq_prof(nlevmax)
1736 
1737       real play(llm), plev_prof_cas(nlev_cas) 
1738       real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
1739       real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
1740       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
1741       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
1742       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
1743       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
1744       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
1745       real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
1746       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
1747 
1748       real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
1749       real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
1750       real u_mod_cas(llm),v_mod_cas(llm)
1751       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
1752       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
1753       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
1754       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
1755       real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
1756       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
1757 
1758       integer l,k,k1,k2
1759       real frac,frac1,frac2,fact
1760 
1761!       do l = 1, llm
1762!       print *,'debut interp2, play=',l,play(l)
1763!       enddo
1764!      do l = 1, nlev_cas
1765!      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
1766!      enddo
1767
1768       do l = 1, llm
1769
1770        if (play(l).ge.plev_prof_cas(nlev_cas)) then
1771 
1772        mxcalc=l
1773!        print *,'debut interp2, mxcalc=',mxcalc
1774         k1=0
1775         k2=0
1776
1777         if (play(l).le.plev_prof_cas(1)) then
1778
1779         do k = 1, nlev_cas-1
1780          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
1781            k1=k
1782            k2=k+1
1783          endif
1784         enddo
1785
1786         if (k1.eq.0 .or. k2.eq.0) then
1787          write(*,*) 'PB! k1, k2 = ',k1,k2
1788          write(*,*) 'l,play(l) = ',l,play(l)/100
1789         do k = 1, nlev_cas-1
1790          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
1791         enddo
1792         endif
1793
1794         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
1795         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
1796         theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
1797         if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
1798         thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
1799         thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
1800         qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
1801         ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
1802         qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
1803         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
1804         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
1805         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
1806         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
1807         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
1808         omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
1809         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
1810         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
1811         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
1812         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
1813         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
1814         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
1815         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
1816         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
1817         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
1818         dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
1819         hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
1820         vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
1821         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
1822         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
1823         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
1824         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
1825     
1826         else !play>plev_prof_cas(1)
1827
1828         k1=1
1829         k2=2
1830         print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
1831         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1832         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1833         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
1834         theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
1835         if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
1836         thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
1837         thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
1838         qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
1839         ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
1840         qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
1841         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
1842         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
1843         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
1844         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
1845         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
1846         omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
1847         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
1848         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
1849         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
1850         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
1851         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
1852         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
1853         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
1854         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
1855         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
1856         dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
1857         hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
1858         vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
1859         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
1860         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
1861         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
1862         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
1863
1864         endif ! play.le.plev_prof_cas(1)
1865
1866        else ! above max altitude of forcing file
1867 
1868!jyg
1869         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
1870         fact = max(fact,0.)                                           !jyg
1871         fact = exp(-fact)                                             !jyg
1872         t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
1873         theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
1874         thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
1875         thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
1876         qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
1877         ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
1878         qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
1879         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
1880         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
1881         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                     !jyg
1882         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                     !jyg
1883         w_mod_cas(l)= 0.0                                             !jyg
1884         omega_mod_cas(l)= 0.0                                         !jyg
1885         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
1886         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
1887         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
1888         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
1889         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
1890         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
1891         dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
1892         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
1893         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
1894         dth_mod_cas(l)= dth_prof_cas(nlev_cas) 
1895         hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
1896         vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
1897         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
1898         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
1899         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
1900         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
1901 
1902        endif ! play
1903 
1904       enddo ! l
1905
1906          return
1907          end
1908!***************************************************************************** 
1909
1910
1911
1912
Note: See TracBrowser for help on using the repository browser.