source: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_1dutils.f90 @ 5158

Last change on this file since 5158 was 5158, checked in by abarral, 7 weeks ago

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables
  • Property svn:keywords set to Id
File size: 61.1 KB
Line 
1MODULE lmdz_1dutils
2  IMPLICIT NONE; PRIVATE
3  PUBLIC fq_sat, conf_unicol, dyn1deta0, dyn1dredem, &
4          disvert0, advect_vert, advect_va, lstendh, nudge_rht_init, nudge_uv_init, &
5          nudge_rht, nudge_uv, interp2_case_vertical
6
7CONTAINS
8  REAL FUNCTION fq_sat(kelvin, millibar)
9    IMPLICIT NONE
10    !======================================================================
11    ! Autheur(s): Z.X. Li (LMD/CNRS)
12    ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
13    !======================================================================
14    ! Arguments:
15    ! kelvin---input-R: temperature en Kelvin
16    ! millibar--input-R: pression en mb
17
18    ! fq_sat----output-R: vapeur d'eau saturante en kg/kg
19    !======================================================================
20
21    REAL, INTENT(IN) :: kelvin, millibar
22
23    REAL r2es
24    PARAMETER (r2es = 611.14 * 18.0153 / 28.9644)
25    REAL r3les, r3ies, r3es
26    PARAMETER (R3LES = 17.269)
27    PARAMETER (R3IES = 21.875)
28
29    REAL r4les, r4ies, r4es
30    PARAMETER (R4LES = 35.86)
31    PARAMETER (R4IES = 7.66)
32
33    REAL rtt
34    PARAMETER (rtt = 273.16)
35
36    REAL retv
37    PARAMETER (retv = 28.9644 / 18.0153 - 1.0)
38
39    REAL zqsat
40    REAL temp, pres
41    !     ------------------------------------------------------------------
42
43    temp = kelvin
44    pres = millibar * 100.0
45    !      WRITE(*,*)'kelvin,millibar=',kelvin,millibar
46    !      WRITE(*,*)'temp,pres=',temp,pres
47
48    IF (temp <= rtt) THEN
49      r3es = r3ies
50      r4es = r4ies
51    ELSE
52      r3es = r3les
53      r4es = r4les
54    ENDIF
55
56    zqsat = r2es / pres * EXP (r3es * (temp - rtt) / (temp - r4es))
57    zqsat = MIN(0.5, ZQSAT)
58    zqsat = zqsat / (1. - retv * zqsat)
59
60    fq_sat = zqsat
61  END FUNCTION fq_sat
62
63  SUBROUTINE conf_unicol
64
65    USE IOIPSL
66    USE lmdz_print_control, ONLY: lunout
67    USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
68    USE lmdz_fcs_gcssold, ONLY: imp_fcg_gcssold, ts_fcg_gcssold, Tp_fcg_gcssold, Tp_ini_gcssold, xTurb_fcg_gcssold
69    USE lmdz_tsoilnudge, ONLY: nudge_tsoil, isoil_nudge, Tsoil_nudge, tau_soil_nudge
70    USE lmdz_compar1d
71
72    IMPLICIT NONE
73    !-----------------------------------------------------------------------
74    !     Auteurs :   A. Lahellec  .
75
76    !  -------------------------------------------------------------------
77
78    !      .........    Initilisation parametres du lmdz1D      ..........
79
80    !---------------------------------------------------------------------
81    !   initialisations:
82    !   ----------------
83
84    !Config  Key  = lunout
85    !Config  Desc = unite de fichier pour les impressions
86    !Config  Def  = 6
87    !Config  Help = unite de fichier pour les impressions
88    !Config         (defaut sortie standard = 6)
89    lunout = 6
90    !      CALL getin('lunout', lunout)
91    IF (lunout /= 5 .AND. lunout /= 6) THEN
92      OPEN(lunout, FILE = 'lmdz.out')
93    ENDIF
94
95    !Config  Key  = prt_level
96    !Config  Desc = niveau d'impressions de debogage
97    !Config  Def  = 0
98    !Config  Help = Niveau d'impression pour le debogage
99    !Config         (0 = minimum d'impression)
100    !      prt_level = 0
101    !      CALL getin('prt_level',prt_level)
102
103    !-----------------------------------------------------------------------
104    !  Parametres de controle du run:
105    !-----------------------------------------------------------------------
106
107    !Config  Key  = restart
108    !Config  Desc = on repart des startphy et start1dyn
109    !Config  Def  = false
110    !Config  Help = les fichiers restart doivent etre renomme en start
111    restart = .FALSE.
112    CALL getin('restart', restart)
113
114    !Config  Key  = forcing_type
115    !Config  Desc = defines the way the SCM is forced:
116    !Config  Def  = 0
117    !!Config  Help = 0 ==> forcing_les = .TRUE.
118    !             initial profiles from file prof.inp.001
119    !             no forcing by LS convergence ;
120    !             surface temperature imposed ;
121    !             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
122    !         = 1 ==> forcing_radconv = .TRUE.
123    !             idem forcing_type = 0, but the imposed radiative cooling
124    !             is set to 0 (hence, if iflag_radia=0 in physiq.def,
125    !             then there is no radiative cooling at all)
126    !         = 2 ==> forcing_toga = .TRUE.
127    !             initial profiles from TOGA-COARE IFA files
128    !             LS convergence and SST imposed from TOGA-COARE IFA files
129    !         = 3 ==> forcing_GCM2SCM = .TRUE.
130    !             initial profiles from the GCM output
131    !             LS convergence imposed from the GCM output
132    !         = 4 ==> forcing_twpi = .TRUE.
133    !             initial profiles from TWPICE nc files
134    !             LS convergence and SST imposed from TWPICE nc files
135    !         = 5 ==> forcing_rico = .TRUE.
136    !             initial profiles from RICO idealized
137    !             LS convergence imposed from  RICO (cst)
138    !         = 6 ==> forcing_amma = .TRUE.
139    !         = 10 ==> forcing_case = .TRUE.
140    !             initial profiles from case.nc file
141    !         = 40 ==> forcing_GCSSold = .TRUE.
142    !             initial profile from GCSS file
143    !             LS convergence imposed from GCSS file
144    !         = 50 ==> forcing_fire = .TRUE.
145    !         = 59 ==> forcing_sandu = .TRUE.
146    !             initial profiles from sanduref file: see prof.inp.001
147    !             SST varying with time and divergence constante: see ifa_sanduref.txt file
148    !             Radiation has to be computed interactively
149    !         = 60 ==> forcing_astex = .TRUE.
150    !             initial profiles from file: see prof.inp.001
151    !             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
152    !             Radiation has to be computed interactively
153    !         = 61 ==> forcing_armcu = .TRUE.
154    !             initial profiles from file: see prof.inp.001
155    !             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
156    !             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
157    !             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
158    !             Radiation to be switched off
159    !         > 100 ==> forcing_case = .TRUE. or forcing_case2 = .TRUE.
160    !             initial profiles from case.nc file
161
162    forcing_type = 0
163    CALL getin('forcing_type', forcing_type)
164    imp_fcg_gcssold = .FALSE.
165    ts_fcg_gcssold = .FALSE.
166    Tp_fcg_gcssold = .FALSE.
167    Tp_ini_gcssold = .FALSE.
168    xTurb_fcg_gcssold = .FALSE.
169    IF (forcing_type ==40) THEN
170      CALL getin('imp_fcg', imp_fcg_gcssold)
171      CALL getin('ts_fcg', ts_fcg_gcssold)
172      CALL getin('tp_fcg', Tp_fcg_gcssold)
173      CALL getin('tp_ini', Tp_ini_gcssold)
174      CALL getin('turb_fcg', xTurb_fcg_gcssold)
175    ENDIF
176
177    !Parametres de forcage
178    !Config  Key  = tend_t
179    !Config  Desc = forcage ou non par advection de T
180    !Config  Def  = false
181    !Config  Help = forcage ou non par advection de T
182    tend_t = 0
183    CALL getin('tend_t', tend_t)
184
185    !Config  Key  = tend_q
186    !Config  Desc = forcage ou non par advection de q
187    !Config  Def  = false
188    !Config  Help = forcage ou non par advection de q
189    tend_q = 0
190    CALL getin('tend_q', tend_q)
191
192    !Config  Key  = tend_u
193    !Config  Desc = forcage ou non par advection de u
194    !Config  Def  = false
195    !Config  Help = forcage ou non par advection de u
196    tend_u = 0
197    CALL getin('tend_u', tend_u)
198
199    !Config  Key  = tend_v
200    !Config  Desc = forcage ou non par advection de v
201    !Config  Def  = false
202    !Config  Help = forcage ou non par advection de v
203    tend_v = 0
204    CALL getin('tend_v', tend_v)
205
206    !Config  Key  = tend_w
207    !Config  Desc = forcage ou non par vitesse verticale
208    !Config  Def  = false
209    !Config  Help = forcage ou non par vitesse verticale
210    tend_w = 0
211    CALL getin('tend_w', tend_w)
212
213    !Config  Key  = tend_rayo
214    !Config  Desc = forcage ou non par dtrad
215    !Config  Def  = false
216    !Config  Help = forcage ou non par dtrad
217    tend_rayo = 0
218    CALL getin('tend_rayo', tend_rayo)
219
220
221    !Config  Key  = nudge_t
222    !Config  Desc = constante de nudging de T
223    !Config  Def  = false
224    !Config  Help = constante de nudging de T
225    nudge_t = 0.
226    CALL getin('nudge_t', nudge_t)
227
228    !Config  Key  = nudge_q
229    !Config  Desc = constante de nudging de q
230    !Config  Def  = false
231    !Config  Help = constante de nudging de q
232    nudge_q = 0.
233    CALL getin('nudge_q', nudge_q)
234
235    !Config  Key  = nudge_u
236    !Config  Desc = constante de nudging de u
237    !Config  Def  = false
238    !Config  Help = constante de nudging de u
239    nudge_u = 0.
240    CALL getin('nudge_u', nudge_u)
241
242    !Config  Key  = nudge_v
243    !Config  Desc = constante de nudging de v
244    !Config  Def  = false
245    !Config  Help = constante de nudging de v
246    nudge_v = 0.
247    CALL getin('nudge_v', nudge_v)
248
249    !Config  Key  = nudge_w
250    !Config  Desc = constante de nudging de w
251    !Config  Def  = false
252    !Config  Help = constante de nudging de w
253    nudge_w = 0.
254    CALL getin('nudge_w', nudge_w)
255
256
257    !Config  Key  = iflag_nudge
258    !Config  Desc = atmospheric nudging ttype (decimal code)
259    !Config  Def  = 0
260    !Config  Help = 0 ==> no nudging
261    !  If digit number n of iflag_nudge is set, then nudging of type n is on
262    !  If digit number n of iflag_nudge is not set, then nudging of type n is off
263    !   (digits are numbered from the right)
264    iflag_nudge = 0
265    CALL getin('iflag_nudge', iflag_nudge)
266
267    !Config  Key  = ok_flux_surf
268    !Config  Desc = forcage ou non par les flux de surface
269    !Config  Def  = false
270    !Config  Help = forcage ou non par les flux de surface
271    ok_flux_surf = .FALSE.
272    CALL getin('ok_flux_surf', ok_flux_surf)
273
274    !Config  Key  = ok_forc_tsurf
275    !Config  Desc = forcage ou non par la Ts
276    !Config  Def  = false
277    !Config  Help = forcage ou non par la Ts
278    ok_forc_tsurf = .FALSE.
279    CALL getin('ok_forc_tsurf', ok_forc_tsurf)
280
281    !Config  Key  = ok_prescr_ust
282    !Config  Desc = ustar impose ou non
283    !Config  Def  = false
284    !Config  Help = ustar impose ou non
285    ok_prescr_ust = .FALSE.
286    CALL getin('ok_prescr_ust', ok_prescr_ust)
287
288
289    !Config  Key  = ok_prescr_beta
290    !Config  Desc = betaevap impose ou non
291    !Config  Def  = false
292    !Config  Help = betaevap impose ou non
293    ok_prescr_beta = .FALSE.
294    CALL getin('ok_prescr_beta', ok_prescr_beta)
295
296    !Config  Key  = ok_old_disvert
297    !Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
298    !Config  Def  = false
299    !Config  Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
300    ok_old_disvert = .FALSE.
301    CALL getin('ok_old_disvert', ok_old_disvert)
302
303    !Config  Key  = time_ini
304    !Config  Desc = meaningless in this  case
305    !Config  Def  = 0.
306    !Config  Help =
307    time_ini = 0.
308    CALL getin('time_ini', time_ini)
309
310    !Config  Key  = rlat et rlon
311    !Config  Desc = latitude et longitude
312    !Config  Def  = 0.0  0.0
313    !Config  Help = fixe la position de la colonne
314    xlat = 0.
315    xlon = 0.
316    CALL getin('rlat', xlat)
317    CALL getin('rlon', xlon)
318
319    !Config  Key  = airephy
320    !Config  Desc = Grid cell area
321    !Config  Def  = 1.e11
322    !Config  Help =
323    airefi = 1.e11
324    CALL getin('airephy', airefi)
325
326    !Config  Key  = nat_surf
327    !Config  Desc = surface type
328    !Config  Def  = 0 (ocean)
329    !Config  Help = 0=ocean,1=land,2=glacier,3=banquise
330    nat_surf = 0.
331    CALL getin('nat_surf', nat_surf)
332
333    !Config  Key  = tsurf
334    !Config  Desc = surface temperature
335    !Config  Def  = 290.
336    !Config  Help = surface temperature
337    tsurf = 290.
338    CALL getin('tsurf', tsurf)
339
340    !Config  Key  = psurf
341    !Config  Desc = surface pressure
342    !Config  Def  = 102400.
343    !Config  Help =
344    psurf = 102400.
345    CALL getin('psurf', psurf)
346
347    !Config  Key  = zsurf
348    !Config  Desc = surface altitude
349    !Config  Def  = 0.
350    !Config  Help =
351    zsurf = 0.
352    CALL getin('zsurf', zsurf)
353    ! EV pour accord avec format standard
354    CALL getin('zorog', zsurf)
355
356
357    !Config  Key  = rugos
358    !Config  Desc = coefficient de frottement
359    !Config  Def  = 0.0001
360    !Config  Help = calcul du Cdrag
361    rugos = 0.0001
362    CALL getin('rugos', rugos)
363    ! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0
364    CALL getin('z0', rugos)
365
366    !Config  Key  = rugosh
367    !Config  Desc = coefficient de frottement
368    !Config  Def  = rugos
369    !Config  Help = calcul du Cdrag
370    rugosh = rugos
371    CALL getin('rugosh', rugosh)
372
373
374
375    !Config  Key  = snowmass
376    !Config  Desc = mass de neige de la surface en kg/m2
377    !Config  Def  = 0.0000
378    !Config  Help = snowmass
379    snowmass = 0.0000
380    CALL getin('snowmass', snowmass)
381
382    !Config  Key  = wtsurf et wqsurf
383    !Config  Desc = ???
384    !Config  Def  = 0.0 0.0
385    !Config  Help =
386    wtsurf = 0.0
387    wqsurf = 0.0
388    CALL getin('wtsurf', wtsurf)
389    CALL getin('wqsurf', wqsurf)
390
391    !Config  Key  = albedo
392    !Config  Desc = albedo
393    !Config  Def  = 0.09
394    !Config  Help =
395    albedo = 0.09
396    CALL getin('albedo', albedo)
397
398    !Config  Key  = agesno
399    !Config  Desc = age de la neige
400    !Config  Def  = 30.0
401    !Config  Help =
402    xagesno = 30.0
403    CALL getin('agesno', xagesno)
404
405    !Config  Key  = restart_runoff
406    !Config  Desc = age de la neige
407    !Config  Def  = 30.0
408    !Config  Help =
409    restart_runoff = 0.0
410    CALL getin('restart_runoff', restart_runoff)
411
412    !Config  Key  = qsolinp
413    !Config  Desc = initial bucket water content (kg/m2) when land (5std)
414    !Config  Def  = 30.0
415    !Config  Help =
416    qsolinp = 1.
417    CALL getin('qsolinp', qsolinp)
418
419
420
421    !Config  Key  = betaevap
422    !Config  Desc = beta for actual evaporation when prescribed
423    !Config  Def  = 1.0
424    !Config  Help =
425    betaevap = 1.
426    CALL getin('betaevap', betaevap)
427
428    !Config  Key  = zpicinp
429    !Config  Desc = denivellation orographie
430    !Config  Def  = 0.
431    !Config  Help =  input brise
432    zpicinp = 0.
433    CALL getin('zpicinp', zpicinp)
434    !Config key = nudge_tsoil
435    !Config  Desc = activation of soil temperature nudging
436    !Config  Def  = .FALSE.
437    !Config  Help = ...
438
439    nudge_tsoil = .FALSE.
440    CALL getin('nudge_tsoil', nudge_tsoil)
441
442    !Config key = isoil_nudge
443    !Config  Desc = level number where soil temperature is nudged
444    !Config  Def  = 3
445    !Config  Help = ...
446
447    isoil_nudge = 3
448    CALL getin('isoil_nudge', isoil_nudge)
449
450    !Config key = Tsoil_nudge
451    !Config  Desc = target temperature for tsoil(isoil_nudge)
452    !Config  Def  = 300.
453    !Config  Help = ...
454
455    Tsoil_nudge = 300.
456    CALL getin('Tsoil_nudge', Tsoil_nudge)
457
458    !Config key = tau_soil_nudge
459    !Config  Desc = nudging relaxation time for tsoil
460    !Config  Def  = 3600.
461    !Config  Help = ...
462
463    tau_soil_nudge = 3600.
464    CALL getin('tau_soil_nudge', tau_soil_nudge)
465
466    !----------------------------------------------------------
467    ! Param??tres de for??age pour les forcages communs:
468    ! Pour les forcages communs: ces entiers valent 0 ou 1
469    ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
470    ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale
471    ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv)
472    ! forcages en omega, w, vent geostrophique ou ustar
473    ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging
474    !----------------------------------------------------------
475
476    !Config  Key  = tadv
477    !Config  Desc = forcage ou non par advection totale de T
478    !Config  Def  = false
479    !Config  Help = forcage ou non par advection totale de T
480    tadv = 0
481    CALL getin('tadv', tadv)
482
483    !Config  Key  = tadvv
484    !Config  Desc = forcage ou non par advection verticale de T
485    !Config  Def  = false
486    !Config  Help = forcage ou non par advection verticale de T
487    tadvv = 0
488    CALL getin('tadvv', tadvv)
489
490    !Config  Key  = tadvh
491    !Config  Desc = forcage ou non par advection horizontale de T
492    !Config  Def  = false
493    !Config  Help = forcage ou non par advection horizontale de T
494    tadvh = 0
495    CALL getin('tadvh', tadvh)
496
497    !Config  Key  = thadv
498    !Config  Desc = forcage ou non par advection totale de Theta
499    !Config  Def  = false
500    !Config  Help = forcage ou non par advection totale de Theta
501    thadv = 0
502    CALL getin('thadv', thadv)
503
504    !Config  Key  = thadvv
505    !Config  Desc = forcage ou non par advection verticale de Theta
506    !Config  Def  = false
507    !Config  Help = forcage ou non par advection verticale de Theta
508    thadvv = 0
509    CALL getin('thadvv', thadvv)
510
511    !Config  Key  = thadvh
512    !Config  Desc = forcage ou non par advection horizontale de Theta
513    !Config  Def  = false
514    !Config  Help = forcage ou non par advection horizontale de Theta
515    thadvh = 0
516    CALL getin('thadvh', thadvh)
517
518    !Config  Key  = qadv
519    !Config  Desc = forcage ou non par advection totale de Q
520    !Config  Def  = false
521    !Config  Help = forcage ou non par advection totale de Q
522    qadv = 0
523    CALL getin('qadv', qadv)
524
525    !Config  Key  = qadvv
526    !Config  Desc = forcage ou non par advection verticale de Q
527    !Config  Def  = false
528    !Config  Help = forcage ou non par advection verticale de Q
529    qadvv = 0
530    CALL getin('qadvv', qadvv)
531
532    !Config  Key  = qadvh
533    !Config  Desc = forcage ou non par advection horizontale de Q
534    !Config  Def  = false
535    !Config  Help = forcage ou non par advection horizontale de Q
536    qadvh = 0
537    CALL getin('qadvh', qadvh)
538
539    !Config  Key  = trad
540    !Config  Desc = forcage ou non par tendance radiative
541    !Config  Def  = false
542    !Config  Help = forcage ou non par tendance radiative
543    trad = 0
544    CALL getin('trad', trad)
545
546    !Config  Key  = forc_omega
547    !Config  Desc = forcage ou non par omega
548    !Config  Def  = false
549    !Config  Help = forcage ou non par omega
550    forc_omega = 0
551    CALL getin('forc_omega', forc_omega)
552
553    !Config  Key  = forc_u
554    !Config  Desc = forcage ou non par u
555    !Config  Def  = false
556    !Config  Help = forcage ou non par u
557    forc_u = 0
558    CALL getin('forc_u', forc_u)
559
560    !Config  Key  = forc_v
561    !Config  Desc = forcage ou non par v
562    !Config  Def  = false
563    !Config  Help = forcage ou non par v
564    forc_v = 0
565    CALL getin('forc_v', forc_v)
566    !Config  Key  = forc_w
567    !Config  Desc = forcage ou non par w
568    !Config  Def  = false
569    !Config  Help = forcage ou non par w
570    forc_w = 0
571    CALL getin('forc_w', forc_w)
572
573    !Config  Key  = forc_geo
574    !Config  Desc = forcage ou non par geo
575    !Config  Def  = false
576    !Config  Help = forcage ou non par geo
577    forc_geo = 0
578    CALL getin('forc_geo', forc_geo)
579
580    ! Meme chose que ok_precr_ust
581    !Config  Key  = forc_ustar
582    !Config  Desc = forcage ou non par ustar
583    !Config  Def  = false
584    !Config  Help = forcage ou non par ustar
585    forc_ustar = 0
586    CALL getin('forc_ustar', forc_ustar)
587    IF (forc_ustar == 1) ok_prescr_ust = .TRUE.
588
589
590    !Config  Key  = nudging_u
591    !Config  Desc = forcage ou non par nudging sur u
592    !Config  Def  = false
593    !Config  Help = forcage ou non par nudging sur u
594    nudging_u = 0
595    CALL getin('nudging_u', nudging_u)
596
597    !Config  Key  = nudging_v
598    !Config  Desc = forcage ou non par nudging sur v
599    !Config  Def  = false
600    !Config  Help = forcage ou non par nudging sur v
601    nudging_v = 0
602    CALL getin('nudging_v', nudging_v)
603
604    !Config  Key  = nudging_w
605    !Config  Desc = forcage ou non par nudging sur w
606    !Config  Def  = false
607    !Config  Help = forcage ou non par nudging sur w
608    nudging_w = 0
609    CALL getin('nudging_w', nudging_w)
610
611    ! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT
612    !Config  Key  = nudging_q
613    !Config  Desc = forcage ou non par nudging sur q
614    !Config  Def  = false
615    !Config  Help = forcage ou non par nudging sur q
616    nudging_qv = 0
617    CALL getin('nudging_q', nudging_qv)
618    CALL getin('nudging_qv', nudging_qv)
619
620    p_nudging_u = 11000.
621    p_nudging_v = 11000.
622    p_nudging_t = 11000.
623    p_nudging_qv = 11000.
624    CALL getin('p_nudging_u', p_nudging_u)
625    CALL getin('p_nudging_v', p_nudging_v)
626    CALL getin('p_nudging_t', p_nudging_t)
627    CALL getin('p_nudging_qv', p_nudging_qv)
628
629    !Config  Key  = nudging_t
630    !Config  Desc = forcage ou non par nudging sur t
631    !Config  Def  = false
632    !Config  Help = forcage ou non par nudging sur t
633    nudging_t = 0
634    CALL getin('nudging_t', nudging_t)
635
636    WRITE(lunout, *)' +++++++++++++++++++++++++++++++++++++++'
637    WRITE(lunout, *)' Configuration des parametres du gcm1D: '
638    WRITE(lunout, *)' +++++++++++++++++++++++++++++++++++++++'
639    WRITE(lunout, *)' restart = ', restart
640    WRITE(lunout, *)' forcing_type = ', forcing_type
641    WRITE(lunout, *)' time_ini = ', time_ini
642    WRITE(lunout, *)' rlat = ', xlat
643    WRITE(lunout, *)' rlon = ', xlon
644    WRITE(lunout, *)' airephy = ', airefi
645    WRITE(lunout, *)' nat_surf = ', nat_surf
646    WRITE(lunout, *)' tsurf = ', tsurf
647    WRITE(lunout, *)' psurf = ', psurf
648    WRITE(lunout, *)' zsurf = ', zsurf
649    WRITE(lunout, *)' rugos = ', rugos
650    WRITE(lunout, *)' snowmass=', snowmass
651    WRITE(lunout, *)' wtsurf = ', wtsurf
652    WRITE(lunout, *)' wqsurf = ', wqsurf
653    WRITE(lunout, *)' albedo = ', albedo
654    WRITE(lunout, *)' xagesno = ', xagesno
655    WRITE(lunout, *)' restart_runoff = ', restart_runoff
656    WRITE(lunout, *)' qsolinp = ', qsolinp
657    WRITE(lunout, *)' zpicinp = ', zpicinp
658    WRITE(lunout, *)' nudge_tsoil = ', nudge_tsoil
659    WRITE(lunout, *)' isoil_nudge = ', isoil_nudge
660    WRITE(lunout, *)' Tsoil_nudge = ', Tsoil_nudge
661    WRITE(lunout, *)' tau_soil_nudge = ', tau_soil_nudge
662    WRITE(lunout, *)' tadv =      ', tadv
663    WRITE(lunout, *)' tadvv =     ', tadvv
664    WRITE(lunout, *)' tadvh =     ', tadvh
665    WRITE(lunout, *)' thadv =     ', thadv
666    WRITE(lunout, *)' thadvv =    ', thadvv
667    WRITE(lunout, *)' thadvh =    ', thadvh
668    WRITE(lunout, *)' qadv =      ', qadv
669    WRITE(lunout, *)' qadvv =     ', qadvv
670    WRITE(lunout, *)' qadvh =     ', qadvh
671    WRITE(lunout, *)' trad =      ', trad
672    WRITE(lunout, *)' forc_omega = ', forc_omega
673    WRITE(lunout, *)' forc_w     = ', forc_w
674    WRITE(lunout, *)' forc_geo   = ', forc_geo
675    WRITE(lunout, *)' forc_ustar = ', forc_ustar
676    WRITE(lunout, *)' nudging_u  = ', nudging_u
677    WRITE(lunout, *)' nudging_v  = ', nudging_v
678    WRITE(lunout, *)' nudging_t  = ', nudging_t
679    WRITE(lunout, *)' nudging_qv  = ', nudging_qv
680    IF (forcing_type ==40) THEN
681      WRITE(lunout, *) '--- Forcing type GCSS Old --- with:'
682      WRITE(lunout, *)'imp_fcg', imp_fcg_gcssold
683      WRITE(lunout, *)'ts_fcg', ts_fcg_gcssold
684      WRITE(lunout, *)'tp_fcg', Tp_fcg_gcssold
685      WRITE(lunout, *)'tp_ini', Tp_ini_gcssold
686      WRITE(lunout, *)'xturb_fcg', xTurb_fcg_gcssold
687    ENDIF
688
689    WRITE(lunout, *)' +++++++++++++++++++++++++++++++++++++++'
690    WRITE(lunout, *)
691
692  END SUBROUTINE conf_unicol
693
694
695  SUBROUTINE dyn1deta0(fichnom, plev, play, phi, phis, presnivs, &
696          &                          ucov, vcov, temp, q, omega2)
697    USE dimphy
698    USE lmdz_grid_phy
699    USE lmdz_phys_para
700    USE iophy
701    USE phys_state_var_mod
702    USE iostart
703    USE lmdz_writefield_phy
704    USE infotrac
705    USE control_mod
706    USE comconst_mod, ONLY: im, jm, lllm
707    USE logic_mod, ONLY: fxyhypb, ysinus
708    USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn
709    USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr
710
711    IMPLICIT NONE
712    !=======================================================
713    ! Ecriture du fichier de redemarrage sous format NetCDF
714    !=======================================================
715    !   Declarations:
716    !   -------------
717    include "dimensions.h"
718
719    !   Arguments:
720    !   ----------
721    CHARACTER*(*) fichnom
722    !Al1 plev tronque pour .nc mais plev(klev+1):=0
723    REAL :: plev(klon, klev + 1), play (klon, klev), phi(klon, klev)
724    REAL :: presnivs(klon, klev)
725    REAL :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev)
726    REAL :: q(klon, klev, nqtot), omega2(klon, klev)
727    !      REAL :: ug(klev),vg(klev),fcoriolis
728    REAL :: phis(klon)
729
730    !   Variables locales pour NetCDF:
731    !   ------------------------------
732    INTEGER iq
733    INTEGER length
734    PARAMETER (length = 100)
735    REAL tab_cntrl(length) ! tableau des parametres du run
736    CHARACTER*4 nmq(nqtot)
737    CHARACTER*12 modname
738    CHARACTER*80 abort_message
739    LOGICAL found
740
741    modname = 'dyn1deta0 : '
742    !!      nmq(1)="vap"
743    !!      nmq(2)="cond"
744    !!      do iq=3,nqtot
745    !!        WRITE(nmq(iq),'("tra",i1)') iq-2
746    !!      enddo
747    DO iq = 1, nqtot
748      nmq(iq) = trim(tracers(iq)%name)
749    ENDDO
750    PRINT*, 'in dyn1deta0 ', fichnom, klon, klev, nqtot
751    CALL open_startphy(fichnom)
752    PRINT*, 'after open startphy ', fichnom, nmq
753
754    ! Lecture des parametres de controle:
755    CALL get_var("controle", tab_cntrl)
756
757    im = tab_cntrl(1)
758    jm = tab_cntrl(2)
759    lllm = tab_cntrl(3)
760    day_ref = tab_cntrl(4)
761    annee_ref = tab_cntrl(5)
762    !      rad        = tab_cntrl(6)
763    !      omeg       = tab_cntrl(7)
764    !      g          = tab_cntrl(8)
765    !      cpp        = tab_cntrl(9)
766    !      kappa      = tab_cntrl(10)
767    !      daysec     = tab_cntrl(11)
768    !      dtvr       = tab_cntrl(12)
769    !      etot0      = tab_cntrl(13)
770    !      ptot0      = tab_cntrl(14)
771    !      ztot0      = tab_cntrl(15)
772    !      stot0      = tab_cntrl(16)
773    !      ang0       = tab_cntrl(17)
774    !      pa         = tab_cntrl(18)
775    !      preff      = tab_cntrl(19)
776
777    !      clon       = tab_cntrl(20)
778    !      clat       = tab_cntrl(21)
779    !      grossismx  = tab_cntrl(22)
780    !      grossismy  = tab_cntrl(23)
781
782    IF (tab_cntrl(24)==1.)  THEN
783      fxyhypb = .TRUE.
784      !        dzoomx   = tab_cntrl(25)
785      !        dzoomy   = tab_cntrl(26)
786      !        taux     = tab_cntrl(28)
787      !        tauy     = tab_cntrl(29)
788    ELSE
789      fxyhypb = .FALSE.
790      ysinus = .FALSE.
791      IF(tab_cntrl(27)==1.) ysinus = .TRUE.
792    ENDIF
793
794    day_ini = tab_cntrl(30)
795    itau_dyn = tab_cntrl(31)
796
797    Print*, 'day_ref,annee_ref,day_ini,itau_dyn', day_ref, annee_ref, day_ini, itau_dyn
798
799    !  Lecture des champs
800    CALL get_field("play", play, found)
801    IF (.NOT. found) PRINT*, modname // 'Le champ <Play> est absent'
802    CALL get_field("phi", phi, found)
803    IF (.NOT. found) PRINT*, modname // 'Le champ <Phi> est absent'
804    CALL get_field("phis", phis, found)
805    IF (.NOT. found) PRINT*, modname // 'Le champ <Phis> est absent'
806    CALL get_field("presnivs", presnivs, found)
807    IF (.NOT. found) PRINT*, modname // 'Le champ <Presnivs> est absent'
808    CALL get_field("ucov", ucov, found)
809    IF (.NOT. found) PRINT*, modname // 'Le champ <ucov> est absent'
810    CALL get_field("vcov", vcov, found)
811    IF (.NOT. found) PRINT*, modname // 'Le champ <vcov> est absent'
812    CALL get_field("temp", temp, found)
813    IF (.NOT. found) PRINT*, modname // 'Le champ <temp> est absent'
814    CALL get_field("omega2", omega2, found)
815    IF (.NOT. found) PRINT*, modname // 'Le champ <omega2> est absent'
816    plev(1, klev + 1) = 0.
817    CALL get_field("plev", plev(:, 1:klev), found)
818    IF (.NOT. found) PRINT*, modname // 'Le champ <Plev> est absent'
819
820    Do iq = 1, nqtot
821      CALL get_field("q" // nmq(iq), q(:, :, iq), found)
822      IF (.NOT.found)PRINT*, modname // 'Le champ <q' // nmq // '> est absent'
823    EndDo
824
825    CALL close_startphy
826    PRINT*, ' close startphy', fichnom, play(1, 1), play(1, klev), temp(1, klev)
827  END SUBROUTINE dyn1deta0
828
829
830  SUBROUTINE dyn1dredem(fichnom, plev, play, phi, phis, presnivs, &
831          &                          ucov, vcov, temp, q, omega2)
832    USE dimphy
833    USE lmdz_grid_phy
834    USE lmdz_phys_para
835    USE phys_state_var_mod
836    USE iostart
837    USE infotrac
838    USE control_mod
839    USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
840    USE logic_mod, ONLY: fxyhypb, ysinus
841    USE temps_mod, ONLY: annee_ref, day_end, day_ref, itau_dyn, itaufin
842    USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr
843
844    IMPLICIT NONE
845    !=======================================================
846    ! Ecriture du fichier de redemarrage sous format NetCDF
847    !=======================================================
848    !   Declarations:
849    !   -------------
850    include "dimensions.h"
851
852    !   Arguments:
853    !   ----------
854    CHARACTER*(*) fichnom
855    !Al1 plev tronque pour .nc mais plev(klev+1):=0
856    REAL :: plev(klon, klev), play (klon, klev), phi(klon, klev)
857    REAL :: presnivs(klon, klev)
858    REAL :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev)
859    REAL :: q(klon, klev, nqtot)
860    REAL :: omega2(klon, klev), rho(klon, klev + 1)
861    !      REAL :: ug(klev),vg(klev),fcoriolis
862    REAL :: phis(klon)
863
864    !   Variables locales pour NetCDF:
865    !   ------------------------------
866    INTEGER nid
867    INTEGER ierr
868    INTEGER iq, l
869    INTEGER length
870    PARAMETER (length = 100)
871    REAL tab_cntrl(length) ! tableau des parametres du run
872    CHARACTER*4 nmq(nqtot)
873    CHARACTER*20 modname
874    CHARACTER*80 abort_message
875
876    INTEGER pass
877
878    CALL open_restartphy(fichnom)
879    PRINT*, 'redm1 ', fichnom, klon, klev, nqtot
880    !!      nmq(1)="vap"
881    !!      nmq(2)="cond"
882    !!      nmq(3)="tra1"
883    !!      nmq(4)="tra2"
884    DO iq = 1, nqtot
885      nmq(iq) = trim(tracers(iq)%name)
886    ENDDO
887
888    !     modname = 'dyn1dredem'
889    !     ierr = nf90_open(fichnom, nf90_write, nid)
890    !     IF (ierr .NE. nf90_noerr) THEN
891    !        abort_message="Pb. d ouverture "//fichnom
892    !        CALL abort_gcm('Modele 1D',abort_message,1)
893    !     ENDIF
894
895    DO l = 1, length
896      tab_cntrl(l) = 0.
897    ENDDO
898    tab_cntrl(1) = FLOAT(iim)
899    tab_cntrl(2) = FLOAT(jjm)
900    tab_cntrl(3) = FLOAT(llm)
901    tab_cntrl(4) = FLOAT(day_ref)
902    tab_cntrl(5) = FLOAT(annee_ref)
903    tab_cntrl(6) = rad
904    tab_cntrl(7) = omeg
905    tab_cntrl(8) = g
906    tab_cntrl(9) = cpp
907    tab_cntrl(10) = kappa
908    tab_cntrl(11) = daysec
909    tab_cntrl(12) = dtvr
910    !       tab_cntrl(13) = etot0
911    !       tab_cntrl(14) = ptot0
912    !       tab_cntrl(15) = ztot0
913    !       tab_cntrl(16) = stot0
914    !       tab_cntrl(17) = ang0
915    !       tab_cntrl(18) = pa
916    !       tab_cntrl(19) = preff
917
918    !    .....    parametres  pour le zoom      ......
919
920    !       tab_cntrl(20)  = clon
921    !       tab_cntrl(21)  = clat
922    !       tab_cntrl(22)  = grossismx
923    !       tab_cntrl(23)  = grossismy
924
925    IF (fxyhypb)   THEN
926      tab_cntrl(24) = 1.
927      !       tab_cntrl(25) = dzoomx
928      !       tab_cntrl(26) = dzoomy
929      tab_cntrl(27) = 0.
930      !       tab_cntrl(28) = taux
931      !       tab_cntrl(29) = tauy
932    ELSE
933      tab_cntrl(24) = 0.
934      !       tab_cntrl(25) = dzoomx
935      !       tab_cntrl(26) = dzoomy
936      tab_cntrl(27) = 0.
937      tab_cntrl(28) = 0.
938      tab_cntrl(29) = 0.
939      IF(ysinus)  tab_cntrl(27) = 1.
940    ENDIF
941    !Al1 iday_end -> day_end
942    tab_cntrl(30) = FLOAT(day_end)
943    tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
944
945    DO pass = 1, 2
946      CALL put_var(pass, "controle", "Param. de controle Dyn1D", tab_cntrl)
947
948      !  Ecriture/extension de la coordonnee temps
949
950
951      !  Ecriture des champs
952
953      CALL put_field(pass, "plev", "p interfaces sauf la nulle", plev)
954      CALL put_field(pass, "play", "", play)
955      CALL put_field(pass, "phi", "geopotentielle", phi)
956      CALL put_field(pass, "phis", "geopotentiell de surface", phis)
957      CALL put_field(pass, "presnivs", "", presnivs)
958      CALL put_field(pass, "ucov", "", ucov)
959      CALL put_field(pass, "vcov", "", vcov)
960      CALL put_field(pass, "temp", "", temp)
961      CALL put_field(pass, "omega2", "", omega2)
962
963      Do iq = 1, nqtot
964        CALL put_field(pass, "q" // nmq(iq), "eau vap ou condens et traceurs", &
965                &                                                      q(:, :, iq))
966      EndDo
967      IF (pass==1) CALL enddef_restartphy
968      IF (pass==2) CALL close_restartphy
969
970    ENDDO
971
972  END SUBROUTINE dyn1dredem
973
974
975  SUBROUTINE disvert0(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
976
977    !    Ancienne version disvert dont on a modifie nom pour utiliser
978    !    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
979    !    (MPL 18092012)
980
981    !    Auteur :  P. Le Van .
982
983    IMPLICIT NONE
984
985    include "dimensions.h"
986    include "paramet.h"
987
988    !=======================================================================
989
990
991    !    s = sigma ** kappa   :  coordonnee  verticale
992    !    dsig(l)            : epaisseur de la couche l ds la coord.  s
993    !    sig(l)             : sigma a l'interface des couches l et l-1
994    !    ds(l)              : distance entre les couches l et l-1 en coord.s
995
996    !=======================================================================
997
998    REAL pa, preff
999    REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
1000    REAL presnivs(llm)
1001
1002    !   declarations:
1003    !   -------------
1004
1005    REAL sig(llm + 1), dsig(llm)
1006
1007    INTEGER l
1008    REAL snorm
1009    REAL alpha, beta, gama, delta, deltaz, h
1010    INTEGER np, ierr
1011    REAL pi, x
1012
1013    !-----------------------------------------------------------------------
1014
1015    pi = 2. * ASIN(1.)
1016
1017    OPEN(99, file = 'sigma.def', status = 'old', form = 'formatted', &
1018            &   iostat = ierr)
1019
1020    !-----------------------------------------------------------------------
1021    !   cas 1 on lit les options dans sigma.def:
1022    !   ----------------------------------------
1023
1024    IF (ierr==0) THEN
1025
1026      PRINT*, 'WARNING!!! on lit les options dans sigma.def'
1027      READ(99, *) deltaz
1028      READ(99, *) h
1029      READ(99, *) beta
1030      READ(99, *) gama
1031      READ(99, *) delta
1032      READ(99, *) np
1033      CLOSE(99)
1034      alpha = deltaz / (llm * h)
1035
1036      DO l = 1, llm
1037        dsig(l) = (alpha + (1. - alpha) * exp(-beta * (llm - l))) * &
1038                &          ((tanh(gama * l) / tanh(gama * llm))**np + &
1039                        &            (1. - l / FLOAT(llm)) * delta)
1040      END DO
1041
1042      sig(1) = 1.
1043      DO l = 1, llm - 1
1044        sig(l + 1) = sig(l) * (1. - dsig(l)) / (1. + dsig(l))
1045      END DO
1046      sig(llm + 1) = 0.
1047
1048      DO l = 1, llm
1049        dsig(l) = sig(l) - sig(l + 1)
1050      END DO
1051
1052    ELSE
1053      !-----------------------------------------------------------------------
1054      !   cas 2 ancienne discretisation (LMD5...):
1055      !   ----------------------------------------
1056
1057      PRINT*, 'WARNING!!! Ancienne discretisation verticale'
1058
1059      h = 7.
1060      snorm = 0.
1061      DO l = 1, llm
1062        x = 2. * asin(1.) * (FLOAT(l) - 0.5) / float(llm + 1)
1063        dsig(l) = 1.0 + 7.0 * SIN(x)**2
1064        snorm = snorm + dsig(l)
1065      ENDDO
1066      snorm = 1. / snorm
1067      DO l = 1, llm
1068        dsig(l) = dsig(l) * snorm
1069      ENDDO
1070      sig(llm + 1) = 0.
1071      DO l = llm, 1, -1
1072        sig(l) = sig(l + 1) + dsig(l)
1073      ENDDO
1074
1075    ENDIF
1076
1077    DO l = 1, llm
1078      nivsigs(l) = FLOAT(l)
1079    ENDDO
1080
1081    DO l = 1, llmp1
1082      nivsig(l) = FLOAT(l)
1083    ENDDO
1084
1085    !    ....  Calculs  de ap(l) et de bp(l)  ....
1086    !    .........................................
1087
1088    !   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
1089
1090    bp(llmp1) = 0.
1091
1092    DO l = 1, llm
1093      !c
1094      !cc    ap(l) = 0.
1095      !cc    bp(l) = sig(l)
1096
1097      bp(l) = EXP(1. - 1. / (sig(l) * sig(l)))
1098      ap(l) = pa * (sig(l) - bp(l))
1099
1100    ENDDO
1101    ap(llmp1) = pa * (sig(llmp1) - bp(llmp1))
1102
1103    PRINT *, ' BP '
1104    PRINT *, bp
1105    PRINT *, ' AP '
1106    PRINT *, ap
1107
1108    DO l = 1, llm
1109      dpres(l) = bp(l) - bp(l + 1)
1110      presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l + 1) + bp(l + 1) * preff)
1111    ENDDO
1112
1113    PRINT *, ' PRESNIVS '
1114    PRINT *, presnivs
1115  END SUBROUTINE disvert0
1116
1117  subroutine advect_vert(llm, w, dt, q, plev)
1118    !===============================================================
1119    !   Schema amont pour l'advection verticale en 1D
1120    !   w est la vitesse verticale dp/dt en Pa/s
1121    !   Traitement en volumes finis
1122    !   d / dt ( zm q ) = delta_z ( omega q )
1123    !   d / dt ( zm ) = delta_z ( omega )
1124    !   avec zm = delta_z ( p )
1125    !   si * designe la valeur au pas de temps t+dt
1126    !   zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1127    !   zm*(l) -zm(l) = w(l+1) - w(l)
1128    !   avec w=omega * dt
1129    !---------------------------------------------------------------
1130    IMPLICIT NONE
1131    ! arguments
1132    INTEGER llm
1133    REAL w(llm + 1), q(llm), plev(llm + 1), dt
1134
1135    ! local
1136    INTEGER l
1137    REAL zwq(llm + 1), zm(llm + 1), zw(llm + 1)
1138    REAL qold
1139
1140    !---------------------------------------------------------------
1141
1142    DO l = 1, llm
1143      zw(l) = dt * w(l)
1144      zm(l) = plev(l) - plev(l + 1)
1145      zwq(l) = q(l) * zw(l)
1146    enddo
1147    zwq(llm + 1) = 0.
1148    zw(llm + 1) = 0.
1149
1150    DO l = 1, llm
1151      qold = q(l)
1152      q(l) = (q(l) * zm(l) + zwq(l + 1) - zwq(l)) / (zm(l) + zw(l + 1) - zw(l))
1153      PRINT*, 'ADV Q ', zm(l), zw(l), zwq(l), qold, q(l)
1154    enddo
1155  end SUBROUTINE advect_vert
1156
1157  SUBROUTINE advect_va(llm, omega, d_t_va, d_q_va, d_u_va, d_v_va, &
1158          &                q, temp, u, v, play)
1159    !itlmd
1160    !----------------------------------------------------------------------
1161    !   Calcul de l'advection verticale (ascendance et subsidence) de
1162    !   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1163    !   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1164    !   sans WTG rajouter une advection horizontale
1165    !----------------------------------------------------------------------
1166    USE lmdz_yomcst
1167
1168    IMPLICIT NONE
1169    !        argument
1170    INTEGER llm
1171    real  omega(llm + 1), d_t_va(llm), d_q_va(llm, 3)
1172    real  d_u_va(llm), d_v_va(llm)
1173    real  q(llm, 3), temp(llm)
1174    real  u(llm), v(llm)
1175    real  play(llm)
1176    ! interne
1177    INTEGER l
1178    REAL alpha, omgdown, omgup
1179
1180    DO l = 1, llm
1181      IF(l==1) THEN
1182        !si omgup pour la couche 1, alors tendance nulle
1183        omgdown = max(omega(2), 0.0)
1184        alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1)))
1185        d_t_va(l) = alpha * (omgdown) - omgdown * (temp(l) - temp(l + 1))             &
1186                & / (play(l) - play(l + 1))
1187
1188        d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :)) / (play(l) - play(l + 1))
1189
1190        d_u_va(l) = -omgdown * (u(l) - u(l + 1)) / (play(l) - play(l + 1))
1191        d_v_va(l) = -omgdown * (v(l) - v(l + 1)) / (play(l) - play(l + 1))
1192
1193      elseif(l==llm) THEN
1194        omgup = min(omega(l), 0.0)
1195        alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1)))
1196        d_t_va(l) = alpha * (omgup) - &
1197
1198                !bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1199                &              omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l))
1200        d_q_va(l, :) = -omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l))
1201        d_u_va(l) = -omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l))
1202        d_v_va(l) = -omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l))
1203
1204      else
1205        omgup = min(omega(l), 0.0)
1206        omgdown = max(omega(l + 1), 0.0)
1207        alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1)))
1208        d_t_va(l) = alpha * (omgup + omgdown) - omgdown * (temp(l) - temp(l + 1))       &
1209                & / (play(l) - play(l + 1)) - &
1210                !bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1211                &              omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l))
1212        !      PRINT*, '  ??? '
1213
1214        d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :))                            &
1215                & / (play(l) - play(l + 1)) - &
1216                &              omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l))
1217        d_u_va(l) = -omgdown * (u(l) - u(l + 1))                                  &
1218                & / (play(l) - play(l + 1)) - &
1219                &              omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l))
1220        d_v_va(l) = -omgdown * (v(l) - v(l + 1))                                  &
1221                & / (play(l) - play(l + 1)) - &
1222                &              omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l))
1223
1224      endif
1225
1226    enddo
1227    !fin itlmd
1228  end SUBROUTINE advect_va
1229
1230
1231  SUBROUTINE lstendH(llm, nqtot, omega, d_t_va, d_q_va, q, temp, u, v, play)
1232    !itlmd
1233    !----------------------------------------------------------------------
1234    !   Calcul de l'advection verticale (ascendance et subsidence) de
1235    !   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
1236    !   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
1237    !   sans WTG rajouter une advection horizontale
1238    !----------------------------------------------------------------------
1239    USE lmdz_yomcst
1240
1241    IMPLICIT NONE
1242    !        argument
1243    INTEGER llm, nqtot
1244    real  omega(llm + 1), d_t_va(llm), d_q_va(llm, nqtot)
1245    !        real  d_u_va(llm), d_v_va(llm)
1246    real  q(llm, nqtot), temp(llm)
1247    real  u(llm), v(llm)
1248    real  play(llm)
1249    REAL cor(llm)
1250    !        real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1251    REAL dph(llm), dqdp(llm), dtdp(llm)
1252    ! interne
1253    INTEGER k
1254    REAL omdn, omup
1255
1256    !        dudp=0.
1257    !        dvdp=0.
1258    dqdp = 0.
1259    dtdp = 0.
1260    !        d_u_va=0.
1261    !        d_v_va=0.
1262
1263    cor(:) = rkappa * temp * (1. + q(:, 1) * rv / rd) / (play * (1. + q(:, 1)))
1264
1265    DO k = 2, llm - 1
1266
1267      dph  (k - 1) = (play(k) - play(k - 1))
1268      !       dudp (k-1) = (u   (k  )- u   (k-1  ))/dph(k-1)
1269      !       dvdp (k-1) = (v   (k  )- v   (k-1  ))/dph(k-1)
1270      dqdp (k - 1) = (q   (k, 1) - q   (k - 1, 1)) / dph(k - 1)
1271      dtdp (k - 1) = (temp(k) - temp(k - 1)) / dph(k - 1)
1272
1273    enddo
1274
1275    !      dudp (  llm  ) = dudp ( llm-1 )
1276    !      dvdp (  llm  ) = dvdp ( llm-1 )
1277    dqdp (llm) = dqdp (llm - 1)
1278    dtdp (llm) = dtdp (llm - 1)
1279
1280    DO k = 2, llm - 1
1281      omdn = max(0.0, omega(k + 1))
1282      omup = min(0.0, omega(k))
1283
1284      !      d_u_va(k)  = -omdn*dudp(k)-omup*dudp(k-1)
1285      !      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
1286      d_q_va(k, 1) = -omdn * dqdp(k) - omup * dqdp(k - 1)
1287      d_t_va(k) = -omdn * dtdp(k) - omup * dtdp(k - 1) + (omup + omdn) * cor(k)
1288    enddo
1289
1290    omdn = max(0.0, omega(2))
1291    omup = min(0.0, omega(llm))
1292    !      d_u_va( 1 )   = -omdn*dudp( 1 )
1293    !      d_u_va(llm)   = -omup*dudp(llm)
1294    !      d_v_va( 1 )   = -omdn*dvdp( 1 )
1295    !      d_v_va(llm)   = -omup*dvdp(llm)
1296    d_q_va(1, 1) = -omdn * dqdp(1)
1297    d_q_va(llm, 1) = -omup * dqdp(llm)
1298    d_t_va(1) = -omdn * dtdp(1) + omdn * cor(1)
1299    d_t_va(llm) = -omup * dtdp(llm)!+omup*cor(llm)
1300
1301    !      IF(abs(rlat(1))>10.) THEN
1302    !     Calculate the tendency due agestrophic motions
1303    !      du_age = fcoriolis*(v-vg)
1304    !      dv_age = fcoriolis*(ug-u)
1305    !      endif
1306
1307    !       CALL writefield_phy('d_t_va',d_t_va,llm)
1308  end SUBROUTINE lstendH
1309
1310
1311  SubROUTINE Nudge_RHT_init(paprs, pplay, t, q, t_targ, rh_targ)
1312    ! ========================================================
1313    USE dimphy
1314    USE lmdz_yoethf
1315
1316    USE lmdz_yomcst
1317
1318    IMPLICIT NONE
1319    INCLUDE "FCTTRE.h"
1320
1321    ! ========================================================
1322    REAL paprs(klon, klevp1)
1323    REAL pplay(klon, klev)
1324
1325    !      Variables d'etat
1326    REAL t(klon, klev)
1327    REAL q(klon, klev)
1328
1329    !   Profiles cible
1330    REAL t_targ(klon, klev)
1331    REAL rh_targ(klon, klev)
1332
1333    INTEGER k, i
1334    REAL zx_qs
1335
1336    DO k = 1, klev
1337      DO i = 1, klon
1338        t_targ(i, k) = t(i, k)
1339        IF (t(i, k)<RTT) THEN
1340          zx_qs = qsats(t(i, k)) / (pplay(i, k))
1341        ELSE
1342          zx_qs = qsatl(t(i, k)) / (pplay(i, k))
1343        ENDIF
1344        rh_targ(i, k) = q(i, k) / zx_qs
1345      ENDDO
1346    ENDDO
1347    print *, 't_targ', t_targ
1348    print *, 'rh_targ', rh_targ
1349
1350  END SUBROUTINE nudge_rht_init
1351
1352  SubROUTINE Nudge_UV_init(paprs, pplay, u, v, u_targ, v_targ)
1353    ! ========================================================
1354    USE dimphy
1355
1356    IMPLICIT NONE
1357
1358    ! ========================================================
1359    REAL paprs(klon, klevp1)
1360    REAL pplay(klon, klev)
1361
1362    !      Variables d'etat
1363    REAL u(klon, klev)
1364    REAL v(klon, klev)
1365
1366    !   Profiles cible
1367    REAL u_targ(klon, klev)
1368    REAL v_targ(klon, klev)
1369
1370    INTEGER k, i
1371
1372    DO k = 1, klev
1373      DO i = 1, klon
1374        u_targ(i, k) = u(i, k)
1375        v_targ(i, k) = v(i, k)
1376      ENDDO
1377    ENDDO
1378    print *, 'u_targ', u_targ
1379    print *, 'v_targ', v_targ
1380
1381    RETURN
1382  END
1383
1384  SubROUTINE Nudge_RHT(dtime, paprs, pplay, t_targ, rh_targ, t, q, &
1385          &                      d_t, d_q)
1386    ! ========================================================
1387    USE dimphy
1388    USE lmdz_yoethf
1389
1390    USE lmdz_yomcst
1391
1392    IMPLICIT NONE
1393    INCLUDE "FCTTRE.h"
1394
1395    ! ========================================================
1396    REAL dtime
1397    REAL paprs(klon, klevp1)
1398    REAL pplay(klon, klev)
1399
1400    !      Variables d'etat
1401    REAL t(klon, klev)
1402    REAL q(klon, klev)
1403
1404    ! Tendances
1405    REAL d_t(klon, klev)
1406    REAL d_q(klon, klev)
1407
1408    !   Profiles cible
1409    REAL t_targ(klon, klev)
1410    REAL rh_targ(klon, klev)
1411
1412    !   Temps de relaxation
1413    REAL tau
1414    !c      DATA tau /3600./
1415    !!      DATA tau /5400./
1416    DATA tau /1800./
1417
1418    INTEGER k, i
1419    REAL zx_qs, rh, tnew, d_rh, rhnew
1420
1421    print *, 'dtime, tau ', dtime, tau
1422    print *, 't_targ', t_targ
1423    print *, 'rh_targ', rh_targ
1424    print *, 'temp ', t
1425    print *, 'hum ', q
1426
1427    DO k = 1, klev
1428      DO i = 1, klon
1429        IF (paprs(i, 1) - pplay(i, k) > 10000.) THEN
1430          IF (t(i, k)<RTT) THEN
1431            zx_qs = qsats(t(i, k)) / (pplay(i, k))
1432          ELSE
1433            zx_qs = qsatl(t(i, k)) / (pplay(i, k))
1434          ENDIF
1435          rh = q(i, k) / zx_qs
1436
1437          d_t(i, k) = d_t(i, k) + 1. / tau * (t_targ(i, k) - t(i, k))
1438          d_rh = 1. / tau * (rh_targ(i, k) - rh)
1439
1440          tnew = t(i, k) + d_t(i, k) * dtime
1441          !jyg<
1442          !   Formule pour q :
1443          !                         d_q = (1/tau) [rh_targ*qsat(T_new) - q]
1444
1445          !  Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
1446          !   qui n'etait pas correcte.
1447
1448          IF (tnew<RTT) THEN
1449            zx_qs = qsats(tnew) / (pplay(i, k))
1450          ELSE
1451            zx_qs = qsatl(tnew) / (pplay(i, k))
1452          ENDIF
1453          !!            d_q(i,k) = d_q(i,k) + d_rh*zx_qs
1454          d_q(i, k) = d_q(i, k) + (1. / tau) * (rh_targ(i, k) * zx_qs - q(i, k))
1455          rhnew = (q(i, k) + d_q(i, k) * dtime) / zx_qs
1456
1457          print *, ' k,d_t,rh,d_rh,rhnew,d_q ', &
1458                  k, d_t(i, k), rh, d_rh, rhnew, d_q(i, k)
1459        ENDIF
1460
1461      ENDDO
1462    ENDDO
1463
1464    RETURN
1465  END
1466
1467  SubROUTINE Nudge_UV(dtime, paprs, pplay, u_targ, v_targ, u, v, &
1468          &                      d_u, d_v)
1469    ! ========================================================
1470    USE dimphy
1471
1472    IMPLICIT NONE
1473
1474    ! ========================================================
1475    REAL dtime
1476    REAL paprs(klon, klevp1)
1477    REAL pplay(klon, klev)
1478
1479    !      Variables d'etat
1480    REAL u(klon, klev)
1481    REAL v(klon, klev)
1482
1483    ! Tendances
1484    REAL d_u(klon, klev)
1485    REAL d_v(klon, klev)
1486
1487    !   Profiles cible
1488    REAL u_targ(klon, klev)
1489    REAL v_targ(klon, klev)
1490
1491    !   Temps de relaxation
1492    REAL tau
1493    !c      DATA tau /3600./
1494    !      DATA tau /5400./
1495    DATA tau /43200./
1496
1497    INTEGER k, i
1498
1499    !print *,'dtime, tau ',dtime,tau
1500    !print *, 'u_targ',u_targ
1501    !print *, 'v_targ',v_targ
1502    !print *,'zonal velocity ',u
1503    !print *,'meridional velocity ',v
1504    DO k = 1, klev
1505      DO i = 1, klon
1506        !CR: nudging everywhere
1507        !           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
1508
1509        d_u(i, k) = d_u(i, k) + 1. / tau * (u_targ(i, k) - u(i, k))
1510        d_v(i, k) = d_v(i, k) + 1. / tau * (v_targ(i, k) - v(i, k))
1511
1512        !           print *,' k,u,d_u,v,d_v ',    &
1513        !                     k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
1514        !           ENDIF
1515
1516      ENDDO
1517    ENDDO
1518
1519    RETURN
1520  END
1521
1522  SUBROUTINE interp2_case_vertical(play, nlev_cas, plev_prof_cas                                    &
1523          &, t_prof_cas, th_prof_cas, thv_prof_cas, thl_prof_cas                                       &
1524          &, qv_prof_cas, ql_prof_cas, qi_prof_cas, u_prof_cas, v_prof_cas                              &
1525          &, ug_prof_cas, vg_prof_cas, vitw_prof_cas, omega_prof_cas                                   &
1526          &, du_prof_cas, hu_prof_cas, vu_prof_cas, dv_prof_cas, hv_prof_cas, vv_prof_cas                &
1527          &, dt_prof_cas, ht_prof_cas, vt_prof_cas, dtrad_prof_cas, dq_prof_cas, hq_prof_cas, vq_prof_cas &
1528          &, dth_prof_cas, hth_prof_cas, vth_prof_cas                                                 &
1529
1530          &, t_mod_cas, theta_mod_cas, thv_mod_cas, thl_mod_cas                                        &
1531          &, qv_mod_cas, ql_mod_cas, qi_mod_cas, u_mod_cas, v_mod_cas                                   &
1532          &, ug_mod_cas, vg_mod_cas, w_mod_cas, omega_mod_cas                                          &
1533          &, du_mod_cas, hu_mod_cas, vu_mod_cas, dv_mod_cas, hv_mod_cas, vv_mod_cas                      &
1534          &, dt_mod_cas, ht_mod_cas, vt_mod_cas, dtrad_mod_cas, dq_mod_cas, hq_mod_cas, vq_mod_cas        &
1535          &, dth_mod_cas, hth_mod_cas, vth_mod_cas, mxcalc)
1536
1537    USE lmdz_yomcst
1538
1539    IMPLICIT NONE
1540
1541    include "dimensions.h"
1542
1543    !-------------------------------------------------------------------------
1544    ! Vertical interpolation of generic case forcing data onto mod_casel levels
1545    !-------------------------------------------------------------------------
1546
1547    INTEGER nlevmax
1548    parameter (nlevmax = 41)
1549    INTEGER nlev_cas, mxcalc
1550    !       real play(llm), plev_prof(nlevmax)
1551    !       real t_prof(nlevmax),q_prof(nlevmax)
1552    !       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1553    !       real ht_prof(nlevmax),vt_prof(nlevmax)
1554    !       real hq_prof(nlevmax),vq_prof(nlevmax)
1555
1556    REAL play(llm), plev_prof_cas(nlev_cas)
1557    REAL t_prof_cas(nlev_cas), th_prof_cas(nlev_cas), thv_prof_cas(nlev_cas), thl_prof_cas(nlev_cas)
1558    REAL qv_prof_cas(nlev_cas), ql_prof_cas(nlev_cas), qi_prof_cas(nlev_cas)
1559    REAL u_prof_cas(nlev_cas), v_prof_cas(nlev_cas)
1560    REAL ug_prof_cas(nlev_cas), vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas), omega_prof_cas(nlev_cas)
1561    REAL du_prof_cas(nlev_cas), hu_prof_cas(nlev_cas), vu_prof_cas(nlev_cas)
1562    REAL dv_prof_cas(nlev_cas), hv_prof_cas(nlev_cas), vv_prof_cas(nlev_cas)
1563    REAL dt_prof_cas(nlev_cas), ht_prof_cas(nlev_cas), vt_prof_cas(nlev_cas), dtrad_prof_cas(nlev_cas)
1564    REAL dth_prof_cas(nlev_cas), hth_prof_cas(nlev_cas), vth_prof_cas(nlev_cas)
1565    REAL dq_prof_cas(nlev_cas), hq_prof_cas(nlev_cas), vq_prof_cas(nlev_cas)
1566
1567    REAL t_mod_cas(llm), theta_mod_cas(llm), thv_mod_cas(llm), thl_mod_cas(llm)
1568    REAL qv_mod_cas(llm), ql_mod_cas(llm), qi_mod_cas(llm)
1569    REAL u_mod_cas(llm), v_mod_cas(llm)
1570    REAL ug_mod_cas(llm), vg_mod_cas(llm), w_mod_cas(llm), omega_mod_cas(llm)
1571    REAL du_mod_cas(llm), hu_mod_cas(llm), vu_mod_cas(llm)
1572    REAL dv_mod_cas(llm), hv_mod_cas(llm), vv_mod_cas(llm)
1573    REAL dt_mod_cas(llm), ht_mod_cas(llm), vt_mod_cas(llm), dtrad_mod_cas(llm)
1574    REAL dth_mod_cas(llm), hth_mod_cas(llm), vth_mod_cas(llm)
1575    REAL dq_mod_cas(llm), hq_mod_cas(llm), vq_mod_cas(llm)
1576
1577    INTEGER l, k, k1, k2
1578    REAL frac, frac1, frac2, fact
1579
1580    !       do l = 1, llm
1581    !       print *,'debut interp2, play=',l,play(l)
1582    !       enddo
1583    !      do l = 1, nlev_cas
1584    !      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
1585    !      enddo
1586
1587    DO l = 1, llm
1588
1589      IF (play(l)>=plev_prof_cas(nlev_cas)) THEN
1590        mxcalc = l
1591        !        print *,'debut interp2, mxcalc=',mxcalc
1592        k1 = 0
1593        k2 = 0
1594
1595        IF (play(l)<=plev_prof_cas(1)) THEN
1596          DO k = 1, nlev_cas - 1
1597            IF (play(l)<=plev_prof_cas(k).AND. play(l)>plev_prof_cas(k + 1)) THEN
1598              k1 = k
1599              k2 = k + 1
1600            endif
1601          enddo
1602
1603          IF (k1==0 .OR. k2==0) THEN
1604            WRITE(*, *) 'PB! k1, k2 = ', k1, k2
1605            WRITE(*, *) 'l,play(l) = ', l, play(l) / 100
1606            DO k = 1, nlev_cas - 1
1607              WRITE(*, *) 'k,plev_prof_cas(k) = ', k, plev_prof_cas(k) / 100
1608            enddo
1609          endif
1610
1611          frac = (plev_prof_cas(k2) - play(l)) / (plev_prof_cas(k2) - plev_prof_cas(k1))
1612          t_mod_cas(l) = t_prof_cas(k2) - frac * (t_prof_cas(k2) - t_prof_cas(k1))
1613          theta_mod_cas(l) = th_prof_cas(k2) - frac * (th_prof_cas(k2) - th_prof_cas(k1))
1614          IF(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD)
1615          thv_mod_cas(l) = thv_prof_cas(k2) - frac * (thv_prof_cas(k2) - thv_prof_cas(k1))
1616          thl_mod_cas(l) = thl_prof_cas(k2) - frac * (thl_prof_cas(k2) - thl_prof_cas(k1))
1617          qv_mod_cas(l) = qv_prof_cas(k2) - frac * (qv_prof_cas(k2) - qv_prof_cas(k1))
1618          ql_mod_cas(l) = ql_prof_cas(k2) - frac * (ql_prof_cas(k2) - ql_prof_cas(k1))
1619          qi_mod_cas(l) = qi_prof_cas(k2) - frac * (qi_prof_cas(k2) - qi_prof_cas(k1))
1620          u_mod_cas(l) = u_prof_cas(k2) - frac * (u_prof_cas(k2) - u_prof_cas(k1))
1621          v_mod_cas(l) = v_prof_cas(k2) - frac * (v_prof_cas(k2) - v_prof_cas(k1))
1622          ug_mod_cas(l) = ug_prof_cas(k2) - frac * (ug_prof_cas(k2) - ug_prof_cas(k1))
1623          vg_mod_cas(l) = vg_prof_cas(k2) - frac * (vg_prof_cas(k2) - vg_prof_cas(k1))
1624          w_mod_cas(l) = vitw_prof_cas(k2) - frac * (vitw_prof_cas(k2) - vitw_prof_cas(k1))
1625          omega_mod_cas(l) = omega_prof_cas(k2) - frac * (omega_prof_cas(k2) - omega_prof_cas(k1))
1626          du_mod_cas(l) = du_prof_cas(k2) - frac * (du_prof_cas(k2) - du_prof_cas(k1))
1627          hu_mod_cas(l) = hu_prof_cas(k2) - frac * (hu_prof_cas(k2) - hu_prof_cas(k1))
1628          vu_mod_cas(l) = vu_prof_cas(k2) - frac * (vu_prof_cas(k2) - vu_prof_cas(k1))
1629          dv_mod_cas(l) = dv_prof_cas(k2) - frac * (dv_prof_cas(k2) - dv_prof_cas(k1))
1630          hv_mod_cas(l) = hv_prof_cas(k2) - frac * (hv_prof_cas(k2) - hv_prof_cas(k1))
1631          vv_mod_cas(l) = vv_prof_cas(k2) - frac * (vv_prof_cas(k2) - vv_prof_cas(k1))
1632          dt_mod_cas(l) = dt_prof_cas(k2) - frac * (dt_prof_cas(k2) - dt_prof_cas(k1))
1633          ht_mod_cas(l) = ht_prof_cas(k2) - frac * (ht_prof_cas(k2) - ht_prof_cas(k1))
1634          vt_mod_cas(l) = vt_prof_cas(k2) - frac * (vt_prof_cas(k2) - vt_prof_cas(k1))
1635          dth_mod_cas(l) = dth_prof_cas(k2) - frac * (dth_prof_cas(k2) - dth_prof_cas(k1))
1636          hth_mod_cas(l) = hth_prof_cas(k2) - frac * (hth_prof_cas(k2) - hth_prof_cas(k1))
1637          vth_mod_cas(l) = vth_prof_cas(k2) - frac * (vth_prof_cas(k2) - vth_prof_cas(k1))
1638          dq_mod_cas(l) = dq_prof_cas(k2) - frac * (dq_prof_cas(k2) - dq_prof_cas(k1))
1639          hq_mod_cas(l) = hq_prof_cas(k2) - frac * (hq_prof_cas(k2) - hq_prof_cas(k1))
1640          vq_mod_cas(l) = vq_prof_cas(k2) - frac * (vq_prof_cas(k2) - vq_prof_cas(k1))
1641          dtrad_mod_cas(l) = dtrad_prof_cas(k2) - frac * (dtrad_prof_cas(k2) - dtrad_prof_cas(k1))
1642
1643        else !play>plev_prof_cas(1)
1644
1645          k1 = 1
1646          k2 = 2
1647          print *, 'interp2_vert, k1,k2=', plev_prof_cas(k1), plev_prof_cas(k2)
1648          frac1 = (play(l) - plev_prof_cas(k2)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
1649          frac2 = (play(l) - plev_prof_cas(k1)) / (plev_prof_cas(k1) - plev_prof_cas(k2))
1650          t_mod_cas(l) = frac1 * t_prof_cas(k1) - frac2 * t_prof_cas(k2)
1651          theta_mod_cas(l) = frac1 * th_prof_cas(k1) - frac2 * th_prof_cas(k2)
1652          IF(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD)
1653          thv_mod_cas(l) = frac1 * thv_prof_cas(k1) - frac2 * thv_prof_cas(k2)
1654          thl_mod_cas(l) = frac1 * thl_prof_cas(k1) - frac2 * thl_prof_cas(k2)
1655          qv_mod_cas(l) = frac1 * qv_prof_cas(k1) - frac2 * qv_prof_cas(k2)
1656          ql_mod_cas(l) = frac1 * ql_prof_cas(k1) - frac2 * ql_prof_cas(k2)
1657          qi_mod_cas(l) = frac1 * qi_prof_cas(k1) - frac2 * qi_prof_cas(k2)
1658          u_mod_cas(l) = frac1 * u_prof_cas(k1) - frac2 * u_prof_cas(k2)
1659          v_mod_cas(l) = frac1 * v_prof_cas(k1) - frac2 * v_prof_cas(k2)
1660          ug_mod_cas(l) = frac1 * ug_prof_cas(k1) - frac2 * ug_prof_cas(k2)
1661          vg_mod_cas(l) = frac1 * vg_prof_cas(k1) - frac2 * vg_prof_cas(k2)
1662          w_mod_cas(l) = frac1 * vitw_prof_cas(k1) - frac2 * vitw_prof_cas(k2)
1663          omega_mod_cas(l) = frac1 * omega_prof_cas(k1) - frac2 * omega_prof_cas(k2)
1664          du_mod_cas(l) = frac1 * du_prof_cas(k1) - frac2 * du_prof_cas(k2)
1665          hu_mod_cas(l) = frac1 * hu_prof_cas(k1) - frac2 * hu_prof_cas(k2)
1666          vu_mod_cas(l) = frac1 * vu_prof_cas(k1) - frac2 * vu_prof_cas(k2)
1667          dv_mod_cas(l) = frac1 * dv_prof_cas(k1) - frac2 * dv_prof_cas(k2)
1668          hv_mod_cas(l) = frac1 * hv_prof_cas(k1) - frac2 * hv_prof_cas(k2)
1669          vv_mod_cas(l) = frac1 * vv_prof_cas(k1) - frac2 * vv_prof_cas(k2)
1670          dt_mod_cas(l) = frac1 * dt_prof_cas(k1) - frac2 * dt_prof_cas(k2)
1671          ht_mod_cas(l) = frac1 * ht_prof_cas(k1) - frac2 * ht_prof_cas(k2)
1672          vt_mod_cas(l) = frac1 * vt_prof_cas(k1) - frac2 * vt_prof_cas(k2)
1673          dth_mod_cas(l) = frac1 * dth_prof_cas(k1) - frac2 * dth_prof_cas(k2)
1674          hth_mod_cas(l) = frac1 * hth_prof_cas(k1) - frac2 * hth_prof_cas(k2)
1675          vth_mod_cas(l) = frac1 * vth_prof_cas(k1) - frac2 * vth_prof_cas(k2)
1676          dq_mod_cas(l) = frac1 * dq_prof_cas(k1) - frac2 * dq_prof_cas(k2)
1677          hq_mod_cas(l) = frac1 * hq_prof_cas(k1) - frac2 * hq_prof_cas(k2)
1678          vq_mod_cas(l) = frac1 * vq_prof_cas(k1) - frac2 * vq_prof_cas(k2)
1679          dtrad_mod_cas(l) = frac1 * dtrad_prof_cas(k1) - frac2 * dtrad_prof_cas(k2)
1680
1681        endif ! play.le.plev_prof_cas(1)
1682
1683      else ! above max altitude of forcing file
1684
1685        !jyg
1686        fact = 20. * (plev_prof_cas(nlev_cas) - play(l)) / plev_prof_cas(nlev_cas) !jyg
1687        fact = max(fact, 0.)                                           !jyg
1688        fact = exp(-fact)                                             !jyg
1689        t_mod_cas(l) = t_prof_cas(nlev_cas)                            !jyg
1690        theta_mod_cas(l) = th_prof_cas(nlev_cas)                       !jyg
1691        thv_mod_cas(l) = thv_prof_cas(nlev_cas)                        !jyg
1692        thl_mod_cas(l) = thl_prof_cas(nlev_cas)                        !jyg
1693        qv_mod_cas(l) = qv_prof_cas(nlev_cas) * fact                     !jyg
1694        ql_mod_cas(l) = ql_prof_cas(nlev_cas) * fact                     !jyg
1695        qi_mod_cas(l) = qi_prof_cas(nlev_cas) * fact                     !jyg
1696        u_mod_cas(l) = u_prof_cas(nlev_cas) * fact                       !jyg
1697        v_mod_cas(l) = v_prof_cas(nlev_cas) * fact                       !jyg
1698        ug_mod_cas(l) = ug_prof_cas(nlev_cas) * fact                     !jyg
1699        vg_mod_cas(l) = vg_prof_cas(nlev_cas) * fact                     !jyg
1700        w_mod_cas(l) = 0.0                                             !jyg
1701        omega_mod_cas(l) = 0.0                                         !jyg
1702        du_mod_cas(l) = du_prof_cas(nlev_cas) * fact
1703        hu_mod_cas(l) = hu_prof_cas(nlev_cas) * fact                     !jyg
1704        vu_mod_cas(l) = vu_prof_cas(nlev_cas) * fact                     !jyg
1705        dv_mod_cas(l) = dv_prof_cas(nlev_cas) * fact
1706        hv_mod_cas(l) = hv_prof_cas(nlev_cas) * fact                     !jyg
1707        vv_mod_cas(l) = vv_prof_cas(nlev_cas) * fact                     !jyg
1708        dt_mod_cas(l) = dt_prof_cas(nlev_cas)
1709        ht_mod_cas(l) = ht_prof_cas(nlev_cas)                          !jyg
1710        vt_mod_cas(l) = vt_prof_cas(nlev_cas)                          !jyg
1711        dth_mod_cas(l) = dth_prof_cas(nlev_cas)
1712        hth_mod_cas(l) = hth_prof_cas(nlev_cas)                        !jyg
1713        vth_mod_cas(l) = vth_prof_cas(nlev_cas)                        !jyg
1714        dq_mod_cas(l) = dq_prof_cas(nlev_cas) * fact
1715        hq_mod_cas(l) = hq_prof_cas(nlev_cas) * fact                     !jyg
1716        vq_mod_cas(l) = vq_prof_cas(nlev_cas) * fact                     !jyg
1717        dtrad_mod_cas(l) = dtrad_prof_cas(nlev_cas) * fact               !jyg
1718
1719      endif ! play
1720
1721    enddo ! l
1722
1723    RETURN
1724  end
1725
1726END MODULE lmdz_1dutils
1727
1728SUBROUTINE gr_fi_dyn(nfield, ngrid, im, jm, pfi, pdyn)
1729  USE lmdz_ssum_scopy, ONLY: scopy
1730
1731  IMPLICIT NONE
1732  !=======================================================================
1733  !   passage d'un champ de la grille scalaire a la grille physique
1734  !=======================================================================
1735
1736  !-----------------------------------------------------------------------
1737  !   declarations:
1738  !   -------------
1739
1740  INTEGER im, jm, ngrid, nfield
1741  REAL pdyn(im, jm, nfield)
1742  REAL pfi(ngrid, nfield)
1743
1744  INTEGER i, j, ifield, ig
1745
1746  !-----------------------------------------------------------------------
1747  !   calcul:
1748  !   -------
1749
1750  DO ifield = 1, nfield
1751    !   traitement des poles
1752    DO i = 1, im
1753      pdyn(i, 1, ifield) = pfi(1, ifield)
1754      pdyn(i, jm, ifield) = pfi(ngrid, ifield)
1755    ENDDO
1756
1757    !   traitement des point normaux
1758    DO j = 2, jm - 1
1759      ig = 2 + (j - 2) * (im - 1)
1760      CALL SCOPY(im - 1, pfi(ig, ifield), 1, pdyn(1, j, ifield), 1)
1761      pdyn(im, j, ifield) = pdyn(1, j, ifield)
1762    ENDDO
1763  ENDDO
1764
1765END SUBROUTINE gr_fi_dyn
1766
1767SUBROUTINE gr_dyn_fi(nfield, im, jm, ngrid, pdyn, pfi)
1768  USE lmdz_ssum_scopy, ONLY: scopy
1769
1770  IMPLICIT NONE
1771  !=======================================================================
1772  !   passage d'un champ de la grille scalaire a la grille physique
1773  !=======================================================================
1774
1775  !-----------------------------------------------------------------------
1776  !   declarations:
1777  !   -------------
1778
1779  INTEGER im, jm, ngrid, nfield
1780  REAL pdyn(im, jm, nfield)
1781  REAL pfi(ngrid, nfield)
1782
1783  INTEGER j, ifield, ig
1784
1785  !-----------------------------------------------------------------------
1786  !   calcul:
1787  !   -------
1788
1789  IF(ngrid/=2 + (jm - 2) * (im - 1).AND.ngrid/=1)                          &
1790          &    STOP 'probleme de dim'
1791  !   traitement des poles
1792  CALL SCOPY(nfield, pdyn, im * jm, pfi, ngrid)
1793  CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(ngrid, 1), ngrid)
1794
1795  !   traitement des point normaux
1796  DO ifield = 1, nfield
1797    DO j = 2, jm - 1
1798      ig = 2 + (j - 2) * (im - 1)
1799      CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1)
1800    ENDDO
1801  ENDDO
1802END SUBROUTINE gr_dyn_fi
Note: See TracBrowser for help on using the repository browser.