Context Navigation

lmdz_1dutils.f90 @ 5225

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

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

Download in other formats:

Original Format