Context Navigation

lmdz_1dutils.f90 @ 5136

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

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

Download in other formats:

Original Format