Context Navigation

lmdz_1dutils.f90 @ 5106

Last change on this file since 5106 was 5106, checked in by abarral, 3 months ago
Turn coefils.h into lmdz_coefils.f90 Put filtreg.F90 inside lmdz_filtreg.F90 Turn mod_filtreg_p.F90 into lmdz_filtreg_p.F90 Delete obsolete parafilt.h* (lint) remove spaces between routine name and args
Property svn:keywords set to `Id`
File size: 62.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format