source: LMDZ6/branches/LMDZ-COSP/libf/phylmd/dyn1d/1DUTILS.h @ 3847

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

Retropedalage sur le 1D.
Annulation de toutes les modifs de la veille.
Frederic

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