Context Navigation

lmdz_1dutils.f90 @ 5139

Last change on this file since 5139 was 5139, checked in by abarral, 8 weeks ago
Put nuage.h, flux_arp.h, compbl.h into modules Move unused phylmd/ini_hist* to obsolete
Property svn:keywords set to `Id`
File size: 61.5 KB

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

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

Download in other formats:

Original Format