source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/dyn1d/1DUTILS.h @ 5441

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

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

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