Context Navigation

lmdz_old_lmdz1d.F90 @ 5447

Last change on this file since 5447 was 5228, checked in by abarral, 4 months ago
Merge r5212 r5213
File size: 55.8 KB

Line
1	! $Id: lmdz1d.F90 3540 2019-06-25 14:50:13Z fairhead $
2
3	MODULE lmdz_old_lmdz1d
4	USE lmdz_old_1dconv, ONLY: copie, get_uvd2, get_uvd, rdgrads, spaces
5	USE lmdz_1dutils, ONLY: interp2_case_vertical, nudge_uv, nudge_rht, lstendh, nudge_uv_init, &
6	nudge_rht_init, disvert0
7	PRIVATE ! -- We'd love to put IMPLICIT NONE; here...
8	PUBLIC old_lmdz1d
9	CONTAINS
10
11	SUBROUTINE old_lmdz1d
12
13	USE ioipsl, ONLY: ju2ymds, ymds2ju, ioconf_calendar, getin
14	USE phys_state_var_mod, ONLY: phys_state_var_init, phys_state_var_end, &
15	clwcon, detr_therm, &
16	qsol, fevap, z0m, z0h, agesno, &
17	du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
18	falb_dir, falb_dif, &
19	ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
20	rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
21	solsw, solswfdiff, t_ancien, q_ancien, u_ancien, v_ancien, &
22	wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
23	wake_deltaq, wake_deltat, wake_s, awake_s, wake_dens, &
24	awake_dens, cv_gen, wake_cstar, &
25	zgam, zmax0, zmea, zpic, zsig, &
26	zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
27	ql_ancien, qs_ancien, qbs_ancien, cf_ancien, rvc_ancien, &
28	prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
29	u10m, v10m, ale_wake, ale_bl_stat
30
31	USE dimphy
32	USE surface_data, ONLY: type_ocean, ok_veget
33	USE pbl_surface_mod, ONLY: ftsoil, pbl_surface_init, pbl_surface_final
34	USE fonte_neige_mod, ONLY: fonte_neige_init, fonte_neige_final
35
36	USE lmdz_infotrac
37	USE control_mod
38	USE indice_sol_mod
39	USE phyaqua_mod
40	USE mod_1D_cases_read2
41	USE mod_1D_amma_read
42	USE lmdz_print_control, ONLY: lunout, prt_level
43	USE iniphysiq_mod, ONLY: iniphysiq
44	USE mod_const_mpi, ONLY: comm_lmdz
45	USE physiq_mod, ONLY: physiq
46	USE comvert_mod, ONLY: presnivs, ap, bp, dpres, nivsig, nivsigs, pa, &
47	preff, aps, bps, pseudoalt, scaleheight
48	USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
49	itau_dyn, itau_phy, start_time, year_len
50	USE phys_cal_mod, ONLY: year_len_phys_cal_mod => year_len
51	USE mod_1D_cases_read, ONLY: interp_case_time ! used in included old_1D_read_forc_cases.h
52	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_OUTPUTPHYSSCM
53	USE lmdz_1dutils, ONLY: fq_sat, conf_unicol, dyn1deta0, dyn1dredem
54	USE lmdz_clesphys
55	USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
56	USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
57	USE lmdz_fcs_gcssold, ONLY: imp_fcg_gcssold, ts_fcg_gcssold, Tp_fcg_gcssold, Tp_ini_gcssold, xTurb_fcg_gcssold
58	USE lmdz_tsoilnudge, ONLY: nudge_tsoil, isoil_nudge, Tsoil_nudge, tau_soil_nudge
59	USE lmdz_yomcst
60	USE lmdz_compar1d
61	USE lmdz_date_cas, ONLY: year_ini_cas, mth_ini_cas, day_deb, heure_ini_cas, pdt_cas, day_ju_ini_cas
62
63	USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
64	USE lmdz_conf_gcm, ONLY: conf_gcm
65
66	IMPLICIT NONE
67
68
69	INCLUDE "dimsoil.h"
70
71	!=====================================================================
72	! DECLARATIONS
73	!=====================================================================
74	!---------------------------------------------------------------------
75	! Arguments d' initialisations de la physique (USER DEFINE)
76	!---------------------------------------------------------------------
77
78	INTEGER, parameter :: ngrid = 1
79	REAL :: zcufi = 1.
80	REAL :: zcvfi = 1.
81
82	!- REAL :: nat_surf
83	!- LOGICAL :: ok_flux_surf
84	!- REAL :: fsens
85	!- REAL :: flat
86	!- REAL :: tsurf
87	!- REAL :: rugos
88	!- REAL :: qsol(1:2)
89	!- REAL :: qsurf
90	!- REAL :: psurf
91	!- REAL :: zsurf
92	!- REAL :: albedo
93	!-
94	!- REAL :: time = 0.
95	!- REAL :: time_ini
96	!- REAL :: xlat
97	!- REAL :: xlon
98	!- REAL :: wtsurf
99	!- REAL :: wqsurf
100	!- REAL :: restart_runoff
101	!- REAL :: xagesno
102	!- REAL :: qsolinp
103	!- REAL :: zpicinp
104	!-
105	REAL :: fnday
106	REAL :: day, daytime
107	REAL :: day1
108	REAL :: heure
109	INTEGER :: jour
110	INTEGER :: mois
111	INTEGER :: an
112
113	!---------------------------------------------------------------------
114	! Declarations related to forcing and initial profiles
115	!---------------------------------------------------------------------
116
117	INTEGER :: kmax = llm
118	INTEGER llm700, nq1, nq2
119	INTEGER, PARAMETER :: nlev_max = 1000, nqmx = 1000
120	REAL timestep, frac
121	REAL height(nlev_max), tttprof(nlev_max), qtprof(nlev_max)
122	real uprof(nlev_max), vprof(nlev_max), e12prof(nlev_max)
123	real ugprof(nlev_max), vgprof(nlev_max), wfls(nlev_max)
124	real dqtdxls(nlev_max), dqtdyls(nlev_max)
125	real dqtdtls(nlev_max), thlpcar(nlev_max)
126	real qprof(nlev_max, nqmx)
127
128	! INTEGER :: forcing_type
129	LOGICAL :: forcing_les = .FALSE.
130	LOGICAL :: forcing_armcu = .FALSE.
131	LOGICAL :: forcing_rico = .FALSE.
132	LOGICAL :: forcing_radconv = .FALSE.
133	LOGICAL :: forcing_toga = .FALSE.
134	LOGICAL :: forcing_twpice = .FALSE.
135	LOGICAL :: forcing_amma = .FALSE.
136	LOGICAL :: forcing_dice = .FALSE.
137	LOGICAL :: forcing_gabls4 = .FALSE.
138
139	LOGICAL :: forcing_GCM2SCM = .FALSE.
140	LOGICAL :: forcing_GCSSold = .FALSE.
141	LOGICAL :: forcing_sandu = .FALSE.
142	LOGICAL :: forcing_astex = .FALSE.
143	LOGICAL :: forcing_fire = .FALSE.
144	LOGICAL :: forcing_case = .FALSE.
145	LOGICAL :: forcing_case2 = .FALSE.
146	LOGICAL :: forcing_SCM = .FALSE.
147	INTEGER :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
148	! (cf read_tsurf1d.F)
149
150	REAL wwww
151	!vertical advection computation
152	! real d_t_z(llm), d_q_z(llm)
153	! real d_t_dyn_z(llm), dq_dyn_z(llm)
154	! real zz(llm)
155	! real zfact
156
157	!flag forcings
158	LOGICAL :: nudge_wind = .TRUE.
159	LOGICAL :: nudge_thermo = .FALSE.
160	LOGICAL :: cptadvw = .TRUE.
161	!=====================================================================
162	! DECLARATIONS FOR EACH CASE
163	!=====================================================================
164
165	INCLUDE "old_1D_decl_cases.h"
166
167	!---------------------------------------------------------------------
168	! Declarations related to nudging
169	!---------------------------------------------------------------------
170	INTEGER :: nudge_max
171	parameter (nudge_max = 9)
172	INTEGER :: inudge_RHT = 1
173	INTEGER :: inudge_UV = 2
174	LOGICAL :: nudge(nudge_max)
175	REAL :: t_targ(llm)
176	REAL :: rh_targ(llm)
177	REAL :: u_targ(llm)
178	REAL :: v_targ(llm)
179
180	!---------------------------------------------------------------------
181	! Declarations related to vertical discretization:
182	!---------------------------------------------------------------------
183	REAL :: pzero = 1.e5
184	REAL :: play (llm), zlay (llm), sig_s(llm), plev(llm + 1)
185	REAL :: playd(llm), zlayd(llm), ap_amma(llm + 1), bp_amma(llm + 1)
186
187	!---------------------------------------------------------------------
188	! Declarations related to variables
189	!---------------------------------------------------------------------
190
191	REAL :: phi(llm)
192	REAL :: teta(llm), tetal(llm), temp(llm), u(llm), v(llm), w(llm)
193	REAL rot(1, llm) ! relative vorticity, in s-1
194	REAL :: rlat_rad(1), rlon_rad(1)
195	REAL :: omega(llm + 1), omega2(llm), rho(llm + 1)
196	REAL :: ug(llm), vg(llm), fcoriolis
197	REAL :: sfdt, cfdt
198	REAL :: du_phys(llm), dv_phys(llm), dt_phys(llm)
199	REAL :: dt_dyn(llm)
200	REAL :: dt_cooling(llm), d_t_adv(llm), d_t_nudge(llm)
201	REAL :: d_u_nudge(llm), d_v_nudge(llm)
202	REAL :: du_adv(llm), dv_adv(llm)
203	REAL :: du_age(llm), dv_age(llm)
204	REAL :: alpha
205	REAL :: ttt
206
207	REAL, ALLOCATABLE, DIMENSION(:, :) :: q
208	REAL, ALLOCATABLE, DIMENSION(:, :) :: dq
209	REAL, ALLOCATABLE, DIMENSION(:, :) :: dq_dyn
210	REAL, ALLOCATABLE, DIMENSION(:, :) :: d_q_adv
211	REAL, ALLOCATABLE, DIMENSION(:, :) :: d_q_nudge
212	! REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
213
214	!---------------------------------------------------------------------
215	! Initialization of surface variables
216	!---------------------------------------------------------------------
217	REAL :: run_off_lic_0(1)
218	REAL :: fder(1), snsrf(1, nbsrf), qsurfsrf(1, nbsrf)
219	REAL :: tsoil(1, nsoilmx, nbsrf)
220	! REAL :: agesno(1,nbsrf)
221
222	!---------------------------------------------------------------------
223	! Call to phyredem
224	!---------------------------------------------------------------------
225	LOGICAL :: ok_writedem = .TRUE.
226	REAL :: sollw_in = 0.
227	REAL :: solsw_in = 0.
228
229	!---------------------------------------------------------------------
230	! Call to physiq
231	!---------------------------------------------------------------------
232	LOGICAL :: firstcall = .TRUE.
233	LOGICAL :: lastcall = .FALSE.
234	REAL :: phis(1) = 0.0
235	REAL :: dpsrf(1)
236
237	!---------------------------------------------------------------------
238	! Initializations of boundary conditions
239	!---------------------------------------------------------------------
240	REAL, ALLOCATABLE :: phy_nat (:) ! 0=ocean libre,1=land,2=glacier,3=banquise
241	REAL, ALLOCATABLE :: phy_alb (:) ! Albedo land only (old value condsurf_jyg=0.3)
242	REAL, ALLOCATABLE :: phy_sst (:) ! SST (will not be used; cf read_tsurf1d.F)
243	REAL, ALLOCATABLE :: phy_bil (:) ! Ne sert que pour les slab_ocean
244	REAL, ALLOCATABLE :: phy_rug (:) ! Longueur rugosite utilisee sur land only
245	REAL, ALLOCATABLE :: phy_ice (:) ! Fraction de glace
246	REAL, ALLOCATABLE :: phy_fter(:) ! Fraction de terre
247	REAL, ALLOCATABLE :: phy_foce(:) ! Fraction de ocean
248	REAL, ALLOCATABLE :: phy_fsic(:) ! Fraction de glace
249	REAL, ALLOCATABLE :: phy_flic(:) ! Fraction de glace
250
251	!---------------------------------------------------------------------
252	! Fichiers et d'autres variables
253	!---------------------------------------------------------------------
254	INTEGER :: k, l, i, it = 1, mxcalc
255	INTEGER :: nsrf
256	INTEGER jcode
257	INTEGER read_climoz
258
259	INTEGER :: it_end ! iteration number of the last call
260	!Al1
261	INTEGER ecrit_slab_oc !1=ecrit,-1=lit,0=no file
262	data ecrit_slab_oc/-1/
263
264	! if flag_inhib_forcing = 0, tendencies of forcing are added
265	! <> 0, tendencies of forcing are not added
266	INTEGER :: flag_inhib_forcing = 0
267
268	!=====================================================================
269	! INITIALIZATIONS
270	!=====================================================================
271	du_phys(:) = 0.
272	dv_phys(:) = 0.
273	dt_phys(:) = 0.
274	dt_dyn(:) = 0.
275	dt_cooling(:) = 0.
276	d_t_adv(:) = 0.
277	d_t_nudge(:) = 0.
278	d_u_nudge(:) = 0.
279	d_v_nudge(:) = 0.
280	du_adv(:) = 0.
281	dv_adv(:) = 0.
282	du_age(:) = 0.
283	dv_age(:) = 0.
284
285	! Initialization of Common turb_forcing
286	dtime_frcg = 0.
287	Turb_fcg_gcssold = .FALSE.
288	hthturb_gcssold = 0.
289	hqturb_gcssold = 0.
290
291
292
293
294	!---------------------------------------------------------------------
295	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
296	!---------------------------------------------------------------------
297	!Al1
298	CALL conf_unicol
299	!Al1 moves this gcssold var from common fcg_gcssold to
300	Turb_fcg_gcssold = xTurb_fcg_gcssold
301	! --------------------------------------------------------------------
302	close(1)
303	!Al1
304	WRITE(, ) 'lmdz1d.def lu => unicol.def'
305
306	! forcing_type defines the way the SCM is forced:
307	!forcing_type = 0 ==> forcing_les = .TRUE.
308	! initial profiles from file prof.inp.001
309	! no forcing by LS convergence ;
310	! surface temperature imposed ;
311	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
312	!forcing_type = 1 ==> forcing_radconv = .TRUE.
313	! idem forcing_type = 0, but the imposed radiative cooling
314	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
315	! then there is no radiative cooling at all)
316	!forcing_type = 2 ==> forcing_toga = .TRUE.
317	! initial profiles from TOGA-COARE IFA files
318	! LS convergence and SST imposed from TOGA-COARE IFA files
319	!forcing_type = 3 ==> forcing_GCM2SCM = .TRUE.
320	! initial profiles from the GCM output
321	! LS convergence imposed from the GCM output
322	!forcing_type = 4 ==> forcing_twpice = .TRUE.
323	! initial profiles from TWP-ICE cdf file
324	! LS convergence, omega and SST imposed from TWP-ICE files
325	!forcing_type = 5 ==> forcing_rico = .TRUE.
326	! initial profiles from RICO files
327	! LS convergence imposed from RICO files
328	!forcing_type = 6 ==> forcing_amma = .TRUE.
329	! initial profiles from AMMA nc file
330	! LS convergence, omega and surface fluxes imposed from AMMA file
331	!forcing_type = 7 ==> forcing_dice = .TRUE.
332	! initial profiles and large scale forcings in dice_driver.nc
333	! Different stages: soil model alone, atm. model alone
334	! then both models coupled
335	!forcing_type = 8 ==> forcing_gabls4 = .TRUE.
336	! initial profiles and large scale forcings in gabls4_driver.nc
337	!forcing_type >= 100 ==> forcing_case = .TRUE.
338	! initial profiles and large scale forcings in cas.nc
339	! LS convergence, omega and SST imposed from CINDY-DYNAMO files
340	! 101=cindynamo
341	! 102=bomex
342	!forcing_type >= 100 ==> forcing_case2 = .TRUE.
343	! temporary flag while all the 1D cases are not whith the same cas.nc forcing file
344	! 103=arm_cu2 ie arm_cu with new forcing format
345	! 104=rico2 ie rico with new forcing format
346	!forcing_type = 40 ==> forcing_GCSSold = .TRUE.
347	! initial profile from GCSS file
348	! LS convergence imposed from GCSS file
349	!forcing_type = 50 ==> forcing_fire = .TRUE.
350	! forcing from fire.nc
351	!forcing_type = 59 ==> forcing_sandu = .TRUE.
352	! initial profiles from sanduref file: see prof.inp.001
353	! SST varying with time and divergence constante: see ifa_sanduref.txt file
354	! Radiation has to be computed interactively
355	!forcing_type = 60 ==> forcing_astex = .TRUE.
356	! initial profiles from file: see prof.inp.001
357	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
358	! Radiation has to be computed interactively
359	!forcing_type = 61 ==> forcing_armcu = .TRUE.
360	! initial profiles from file: see prof.inp.001
361	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
362	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
363	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
364	! Radiation to be switched off
365
366	IF (forcing_type <=0) THEN
367	forcing_les = .TRUE.
368	elseif (forcing_type ==1) THEN
369	forcing_radconv = .TRUE.
370	elseif (forcing_type ==2) THEN
371	forcing_toga = .TRUE.
372	elseif (forcing_type ==3) THEN
373	forcing_GCM2SCM = .TRUE.
374	elseif (forcing_type ==4) THEN
375	forcing_twpice = .TRUE.
376	elseif (forcing_type ==5) THEN
377	forcing_rico = .TRUE.
378	elseif (forcing_type ==6) THEN
379	forcing_amma = .TRUE.
380	elseif (forcing_type ==7) THEN
381	forcing_dice = .TRUE.
382	elseif (forcing_type ==8) THEN
383	forcing_gabls4 = .TRUE.
384	elseif (forcing_type ==101) THEN ! Cindynamo starts 1-10-2011 0h
385	forcing_case = .TRUE.
386	year_ini_cas = 2011
387	mth_ini_cas = 10
388	day_deb = 1
389	heure_ini_cas = 0.
390	pdt_cas = 3 * 3600. ! forcing frequency
391	elseif (forcing_type ==102) THEN ! Bomex starts 24-6-1969 0h
392	forcing_case = .TRUE.
393	year_ini_cas = 1969
394	mth_ini_cas = 6
395	day_deb = 24
396	heure_ini_cas = 0.
397	pdt_cas = 1800. ! forcing frequency
398	elseif (forcing_type ==103) THEN ! Arm_cu starts 21-6-1997 11h30
399	forcing_case2 = .TRUE.
400	year_ini_cas = 1997
401	mth_ini_cas = 6
402	day_deb = 21
403	heure_ini_cas = 11.5
404	pdt_cas = 1800. ! forcing frequency
405	elseif (forcing_type ==104) THEN ! rico starts 16-12-2004 0h
406	forcing_case2 = .TRUE.
407	year_ini_cas = 2004
408	mth_ini_cas = 12
409	day_deb = 16
410	heure_ini_cas = 0.
411	pdt_cas = 1800. ! forcing frequency
412	elseif (forcing_type ==105) THEN ! bomex starts 16-12-2004 0h
413	forcing_case2 = .TRUE.
414	year_ini_cas = 1969
415	mth_ini_cas = 6
416	day_deb = 24
417	heure_ini_cas = 0.
418	pdt_cas = 1800. ! forcing frequency
419	elseif (forcing_type ==106) THEN ! ayotte_24SC starts 6-11-1992 0h
420	forcing_case2 = .TRUE.
421	year_ini_cas = 1992
422	mth_ini_cas = 11
423	day_deb = 6
424	heure_ini_cas = 10.
425	pdt_cas = 86400. ! forcing frequency
426	elseif (forcing_type ==113) THEN ! Arm_cu starts 21-6-1997 11h30
427	forcing_SCM = .TRUE.
428	year_ini_cas = 1997
429	! It is possible that those parameters are run twice.
430	CALL getin('anneeref', year_ini_cas)
431	CALL getin('dayref', day_deb)
432	mth_ini_cas = 1 ! pour le moment on compte depuis le debut de l'annee
433	CALL getin('time_ini', heure_ini_cas)
434	elseif (forcing_type ==40) THEN
435	forcing_GCSSold = .TRUE.
436	elseif (forcing_type ==50) THEN
437	forcing_fire = .TRUE.
438	elseif (forcing_type ==59) THEN
439	forcing_sandu = .TRUE.
440	elseif (forcing_type ==60) THEN
441	forcing_astex = .TRUE.
442	elseif (forcing_type ==61) THEN
443	forcing_armcu = .TRUE.
444	IF(llm/=19.AND.llm/=40) stop 'Erreur nombre de niveaux !!'
445	else
446	write (, ) 'ERROR : unknown forcing_type ', forcing_type
447	stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
448	ENDIF
449	PRINT*, "forcing type=", forcing_type
450
451	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
452	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
453	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
454	! through the common sst_forcing.
455
456	type_ts_forcing = 0
457	IF (forcing_toga.OR.forcing_sandu.OR.forcing_astex .OR. forcing_dice) &
458	type_ts_forcing = 1
459
460	! Initialization of the LOGICAL switch for nudging
461	jcode = iflag_nudge
462	DO i = 1, nudge_max
463	nudge(i) = mod(jcode, 10) >= 1
464	jcode = jcode / 10
465	enddo
466	!---------------------------------------------------------------------
467	! Definition of the run
468	!---------------------------------------------------------------------
469
470	CALL conf_gcm(99, .TRUE.)
471
472	!-----------------------------------------------------------------------
473	allocate(phy_nat (year_len)) ! 0=ocean libre,1=land,2=glacier,3=banquise
474	phy_nat(:) = 0.0
475	allocate(phy_alb (year_len)) ! Albedo land only (old value condsurf_jyg=0.3)
476	allocate(phy_sst (year_len)) ! SST (will not be used; cf read_tsurf1d.F)
477	allocate(phy_bil (year_len)) ! Ne sert que pour les slab_ocean
478	phy_bil(:) = 1.0
479	allocate(phy_rug (year_len)) ! Longueur rugosite utilisee sur land only
480	allocate(phy_ice (year_len)) ! Fraction de glace
481	phy_ice(:) = 0.0
482	allocate(phy_fter(year_len)) ! Fraction de terre
483	phy_fter(:) = 0.0
484	allocate(phy_foce(year_len)) ! Fraction de ocean
485	phy_foce(:) = 0.0
486	allocate(phy_fsic(year_len)) ! Fraction de glace
487	phy_fsic(:) = 0.0
488	allocate(phy_flic(year_len)) ! Fraction de glace
489	phy_flic(:) = 0.0
490	!-----------------------------------------------------------------------
491	! Choix du calendrier
492	! -------------------
493
494	! calend = 'earth_365d'
495	IF (calend == 'earth_360d') THEN
496	CALL ioconf_calendar('360_day')
497	WRITE(, )'CALENDRIER CHOISI: Terrestre a 360 jours/an'
498	ELSE IF (calend == 'earth_365d') THEN
499	CALL ioconf_calendar('noleap')
500	WRITE(, )'CALENDRIER CHOISI: Terrestre a 365 jours/an'
501	ELSE IF (calend == 'earth_366d') THEN
502	CALL ioconf_calendar('all_leap')
503	WRITE(, )'CALENDRIER CHOISI: Terrestre bissextile'
504	ELSE IF (calend == 'gregorian') THEN
505	stop 'gregorian calend should not be used by normal user'
506	CALL ioconf_calendar('gregorian') ! not to be used by normal users
507	WRITE(, )'CALENDRIER CHOISI: Gregorien'
508	else
509	write (, ) 'ERROR : unknown calendar ', calend
510	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
511	endif
512	!-----------------------------------------------------------------------
513
514	!c Date :
515	! La date est supposee donnee sous la forme [annee, numero du jour dans
516	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
517	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
518	! Le numero du jour est dans "day". L heure est traitee separement.
519	! La date complete est dans "daytime" (l'unite est le jour).
520	IF (nday>0) THEN
521	fnday = nday
522	else
523	fnday = -nday / float(day_step)
524	endif
525	PRINT *, 'fnday=', fnday
526	! start_time doit etre en FRACTION DE JOUR
527	start_time = time_ini / 24.
528
529	! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
530	IF(forcing_type == 61) fnday = 53100. / 86400.
531	IF(forcing_type == 103) fnday = 53100. / 86400.
532	! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
533	IF(forcing_type == 6) fnday = 64800. / 86400.
534	! IF(forcing_type .EQ. 6) fnday=50400./86400.
535	IF(forcing_type == 8) fnday = 129600. / 86400.
536	annee_ref = anneeref
537	mois = 1
538	day_ref = dayref
539	heure = 0.
540	itau_dyn = 0
541	itau_phy = 0
542	CALL ymds2ju(annee_ref, mois, day_ref, heure, day)
543	day_ini = int(day)
544	day_end = day_ini + int(fnday)
545
546	IF (forcing_type ==2) THEN
547	! Convert the initial date of Toga-Coare to Julian day
548	CALL ymds2ju &
549	(year_ini_toga, mth_ini_toga, day_ini_toga, heure, day_ju_ini_toga)
550
551	ELSEIF (forcing_type ==4) THEN
552	! Convert the initial date of TWPICE to Julian day
553	CALL ymds2ju &
554	(year_ini_twpi, mth_ini_twpi, day_ini_twpi, heure_ini_twpi &
555	, day_ju_ini_twpi)
556	ELSEIF (forcing_type ==6) THEN
557	! Convert the initial date of AMMA to Julian day
558	CALL ymds2ju &
559	(year_ini_amma, mth_ini_amma, day_ini_amma, heure_ini_amma &
560	, day_ju_ini_amma)
561	ELSEIF (forcing_type ==7) THEN
562	! Convert the initial date of DICE to Julian day
563	CALL ymds2ju &
564	(year_ini_dice, mth_ini_dice, day_ini_dice, heure_ini_dice &
565	, day_ju_ini_dice)
566	ELSEIF (forcing_type ==8) THEN
567	! Convert the initial date of GABLS4 to Julian day
568	CALL ymds2ju &
569	(year_ini_gabls4, mth_ini_gabls4, day_ini_gabls4, heure_ini_gabls4 &
570	, day_ju_ini_gabls4)
571	ELSEIF (forcing_type >100) THEN
572	! Convert the initial date to Julian day
573	day_ini_cas = day_deb
574	PRINT*, 'time case', year_ini_cas, mth_ini_cas, day_ini_cas
575	CALL ymds2ju &
576	(year_ini_cas, mth_ini_cas, day_ini_cas, heure_ini_cas * 3600 &
577	, day_ju_ini_cas)
578	PRINT*, 'time case 2', day_ini_cas, day_ju_ini_cas
579	ELSEIF (forcing_type ==59) THEN
580	! Convert the initial date of Sandu case to Julian day
581	CALL ymds2ju &
582	(year_ini_sandu, mth_ini_sandu, day_ini_sandu, &
583	time_ini * 3600., day_ju_ini_sandu)
584
585	ELSEIF (forcing_type ==60) THEN
586	! Convert the initial date of Astex case to Julian day
587	CALL ymds2ju &
588	(year_ini_astex, mth_ini_astex, day_ini_astex, &
589	time_ini * 3600., day_ju_ini_astex)
590
591	ELSEIF (forcing_type ==61) THEN
592	! Convert the initial date of Arm_cu case to Julian day
593	CALL ymds2ju &
594	(year_ini_armcu, mth_ini_armcu, day_ini_armcu, heure_ini_armcu &
595	, day_ju_ini_armcu)
596	ENDIF
597
598	IF (forcing_type >100) THEN
599	daytime = day + heure_ini_cas / 24. ! 1st day and initial time of the simulation
600	ELSE
601	daytime = day + time_ini / 24. ! 1st day and initial time of the simulation
602	ENDIF
603	! Print out the actual date of the beginning of the simulation :
604	CALL ju2ymds(daytime, year_print, month_print, day_print, sec_print)
605	PRINT *, ' Time of beginning : ', &
606	year_print, month_print, day_print, sec_print
607
608	!---------------------------------------------------------------------
609	! Initialization of dimensions, geometry and initial state
610	!---------------------------------------------------------------------
611	! CALL init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
612	! but we still need to initialize dimphy module (klon,klev,etc.) here.
613	CALL init_dimphy1D(1, llm)
614	CALL suphel
615	CALL init_infotrac
616
617	IF (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
618	allocate(q(llm, nqtot)) ; q(:, :) = 0.
619	allocate(dq(llm, nqtot))
620	allocate(dq_dyn(llm, nqtot))
621	allocate(d_q_adv(llm, nqtot))
622	allocate(d_q_nudge(llm, nqtot))
623	! allocate(d_th_adv(llm))
624
625	q(:, :) = 0.
626	dq(:, :) = 0.
627	dq_dyn(:, :) = 0.
628	d_q_adv(:, :) = 0.
629	d_q_nudge(:, :) = 0.
630
631	! No ozone climatology need be read in this pre-initialization
632	! (phys_state_var_init is called again in physiq)
633	read_climoz = 0
634	nsw = 6 ! EV et LF: sinon, falb_dir et falb_dif ne peuvent etre alloues
635
636	CALL phys_state_var_init(read_climoz)
637
638	IF (ngrid/=klon) THEN
639	PRINT*, 'stop in inifis'
640	PRINT*, 'Probleme de dimensions :'
641	PRINT*, 'ngrid = ', ngrid
642	PRINT*, 'klon = ', klon
643	stop
644	endif
645	!!!=====================================================================
646	!!! Feedback forcing values for Gateaux differentiation (al1)
647	!!!=====================================================================
648	!!! Surface Planck forcing bracketing CALL radiation
649	!! surf_Planck = 0.
650	!! surf_Conv = 0.
651	!! WRITE(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
652	!!! a mettre dans le lmdz1d.def ou autre
653	!!
654	!!
655	qsol = qsolinp
656	qsurf = fq_sat(tsurf, psurf / 100.)
657	beta_surf = 1.
658	beta_aridity(:, :) = beta_surf
659	day1 = day_ini
660	time = daytime - day
661	ts_toga(1) = tsurf ! needed by read_tsurf1d.F
662	rho(1) = psurf / (rd * tsurf * (1. + (rv / rd - 1.) * qsurf))
663
664	!! mpl et jyg le 22/08/2012 :
665	!! pour que les cas a flux de surface imposes marchent
666	IF(.NOT.ok_flux_surf.OR.max(abs(wtsurf), abs(wqsurf))>0.) THEN
667	fsens = -wtsurf * rcpd * rho(1)
668	flat = -wqsurf * rlvtt * rho(1)
669	PRINT *, 'Flux: ok_flux wtsurf wqsurf', ok_flux_surf, wtsurf, wqsurf
670	ENDIF
671	PRINT*, 'Flux sol ', fsens, flat
672	!! ok_flux_surf=.FALSE.
673	!! fsens=-wtsurfrcpdrho(1)
674	!! flat=-wqsurfrlvttrho(1)
675	!!!!
676
677	! Vertical discretization and pressure levels at half and mid levels:
678
679	pa = 5e4
680	!! preff= 1.01325e5
681	preff = psurf
682	IF (ok_old_disvert) THEN
683	CALL disvert0(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
684	PRINT *, 'On utilise disvert0'
685	aps(1:llm) = 0.5 * (ap(1:llm) + ap(2:llm + 1))
686	bps(1:llm) = 0.5 * (bp(1:llm) + bp(2:llm + 1))
687	scaleheight = 8.
688	pseudoalt(1:llm) = -scaleheight * log(presnivs(1:llm) / preff)
689	ELSE
690	CALL disvert()
691	PRINT *, 'On utilise disvert'
692	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
693	! Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
694	ENDIF
695
696	sig_s = presnivs / preff
697	plev = ap + bp * psurf
698	play = 0.5 * (plev(1:llm) + plev(2:llm + 1))
699	zlay = -rd * 300. * log(play / psurf) / rg ! moved after reading profiles
700
701	IF (forcing_type == 59) THEN
702	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
703	WRITE(, ) '***********************'
704	DO l = 1, llm
705	WRITE(, ) 'l,play(l),presnivs(l): ', l, play(l), presnivs(l)
706	IF (trouve_700 .AND. play(l)<=70000) THEN
707	llm700 = l
708	PRINT *, 'llm700,play=', llm700, play(l) / 100.
709	trouve_700 = .FALSE.
710	endif
711	enddo
712	WRITE(, ) '***********************'
713	ENDIF
714
715	!=====================================================================
716	! EVENTUALLY, READ FORCING DATA :
717	!=====================================================================
718
719	INCLUDE "old_1D_read_forc_cases.h"
720
721	IF (forcing_GCM2SCM) THEN
722	write (, ) 'forcing_GCM2SCM not yet implemented'
723	stop 'in initialization'
724	END IF ! forcing_GCM2SCM
725
726	PRINT*, 'mxcalc=', mxcalc
727	! PRINT*,'zlay=',zlay(mxcalc)
728	PRINT*, 'play=', play(mxcalc)
729
730	!Al1 pour SST forced, appell?? depuis ocean_forced_noice
731	! EV tg instead of ts_cur
732
733	tg = tsurf ! SST used in read_tsurf1d
734	!=====================================================================
735	! Initialisation de la physique :
736	!=====================================================================
737
738	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
739
740	! day_step, iphysiq lus dans gcm.def ci-dessus
741	! timestep: calcule ci-dessous from rday et day_step
742	! ngrid=1
743	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
744	! rday: defini dans suphel.F (86400.)
745	! day_ini: lu dans run.def (dayref)
746	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
747	! airefi,zcufi,zcvfi initialises au debut de ce programme
748	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
749	day_step = float(nsplit_phys) * day_step / float(iphysiq)
750	write (, ) 'Time step divided by nsplit_phys (=', nsplit_phys, ')'
751	timestep = rday / day_step
752	dtime_frcg = timestep
753
754	zcufi = airefi
755	zcvfi = airefi
756
757	rlat_rad(1) = xlat * rpi / 180.
758	rlon_rad(1) = xlon * rpi / 180.
759
760	! iniphysiq will CALL iniaqua who needs year_len from phys_cal_mod
761	year_len_phys_cal_mod = year_len
762
763	! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
764	! e.g. for cell boundaries, which are meaningless in 1D; so pad these
765	! with '0.' when necessary
766	CALL iniphysiq(iim, jjm, llm, &
767	1, comm_lmdz, &
768	rday, day_ini, timestep, &
769	(/rlat_rad(1), 0./), (/0./), &
770	(/0., 0./), (/rlon_rad(1), 0./), &
771	(/ (/airefi, 0./), (/0., 0./) /), &
772	(/zcufi, 0., 0., 0./), &
773	(/zcvfi, 0./), &
774	ra, rg, rd, rcpd, 1)
775	PRINT*, 'apres iniphysiq'
776
777	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
778	co2_ppm = 330.0
779	solaire = 1370.0
780
781	! Ecriture du startphy avant le premier appel a la physique.
782	! On le met juste avant pour avoir acces a tous les champs
783
784	IF (ok_writedem) THEN
785	!--------------------------------------------------------------------------
786	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
787	! need : qsol fder snow qsurf evap rugos agesno ftsoil
788	!--------------------------------------------------------------------------
789
790	type_ocean = "force"
791	run_off_lic_0(1) = restart_runoff
792	CALL fonte_neige_init(run_off_lic_0)
793
794	fder = 0.
795	snsrf(1, :) = snowmass ! masse de neige des sous surface
796	PRINT *, 'snsrf', snsrf
797	qsurfsrf(1, :) = qsurf ! humidite de l'air des sous surface
798	fevap = 0.
799	z0m(1, :) = rugos ! couverture de neige des sous surface
800	z0h(1, :) = rugosh ! couverture de neige des sous surface
801	agesno = xagesno
802	tsoil(:, :, :) = tsurf
803	!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
804	! tsoil(1,1,1)=299.18
805	! tsoil(1,2,1)=300.08
806	! tsoil(1,3,1)=301.88
807	! tsoil(1,4,1)=305.48
808	! tsoil(1,5,1)=308.00
809	! tsoil(1,6,1)=308.00
810	! tsoil(1,7,1)=308.00
811	! tsoil(1,8,1)=308.00
812	! tsoil(1,9,1)=308.00
813	! tsoil(1,10,1)=308.00
814	! tsoil(1,11,1)=308.00
815	!-----------------------------------------------------------------------
816	CALL pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
817
818	!------------------ prepare limit conditions for limit.nc -----------------
819	!-- Ocean force
820
821	PRINT*, 'avant phyredem'
822	pctsrf(1, :) = 0.
823	IF (nat_surf==0.) THEN
824	pctsrf(1, is_oce) = 1.
825	pctsrf(1, is_ter) = 0.
826	pctsrf(1, is_lic) = 0.
827	pctsrf(1, is_sic) = 0.
828	ELSE IF (nat_surf == 1) THEN
829	pctsrf(1, is_oce) = 0.
830	pctsrf(1, is_ter) = 1.
831	pctsrf(1, is_lic) = 0.
832	pctsrf(1, is_sic) = 0.
833	ELSE IF (nat_surf == 2) THEN
834	pctsrf(1, is_oce) = 0.
835	pctsrf(1, is_ter) = 0.
836	pctsrf(1, is_lic) = 1.
837	pctsrf(1, is_sic) = 0.
838	ELSE IF (nat_surf == 3) THEN
839	pctsrf(1, is_oce) = 0.
840	pctsrf(1, is_ter) = 0.
841	pctsrf(1, is_lic) = 0.
842	pctsrf(1, is_sic) = 1.
843
844	end if
845
846	PRINT*, 'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)', nat_surf &
847	, pctsrf(1, is_oce), pctsrf(1, is_ter)
848
849	zmasq = pctsrf(1, is_ter) + pctsrf(1, is_lic)
850	zpic = zpicinp
851	ftsol = tsurf
852	falb_dir = albedo
853	falb_dif = albedo
854	rugoro = rugos
855	t_ancien(1, :) = temp(:)
856	q_ancien(1, :) = q(:, 1)
857	ql_ancien = 0.
858	qs_ancien = 0.
859	prlw_ancien = 0.
860	prsw_ancien = 0.
861	prw_ancien = 0.
862	IF ( ok_bs ) THEN
863	qbs_ancien = 0.
864	prbsw_ancien = 0.
865	END IF
866	IF (ok_ice_supersat) THEN
867	cf_ancien = 0.
868	rvc_ancien = 0.
869	END IF
870	!jyg<
871	!! pbl_tke(:,:,:)=1.e-8
872	pbl_tke(:, :, :) = 0.
873	pbl_tke(:, 2, :) = 1.e-2
874	PRINT *, ' pbl_tke dans lmdz1d '
875	IF (prt_level >= 5) THEN
876	DO nsrf = 1, 4
877	PRINT *, 'pbl_tke(1,:,', nsrf, ') ', pbl_tke(1, :, nsrf)
878	ENDDO
879	end if
880
881	!>jyg
882
883	rain_fall = 0.
884	snow_fall = 0.
885	solsw = 0.
886	solswfdiff = 0.
887	sollw = 0.
888	sollwdown = rsigma * tsurf**4
889	radsol = 0.
890	rnebcon = 0.
891	ratqs = 0.
892	clwcon = 0.
893	zmax0 = 0.
894	zmea = zsurf
895	zstd = 0.
896	zsig = 0.
897	zgam = 0.
898	zval = 0.
899	zthe = 0.
900	sig1 = 0.
901	w01 = 0.
902
903	wake_deltaq = 0.
904	wake_deltat = 0.
905	wake_delta_pbl_TKE(:, :, :) = 0.
906	delta_tsurf = 0.
907	wake_fip = 0.
908	wake_pe = 0.
909	wake_s = 0.
910	awake_s = 0.
911	wake_dens = 0.
912	awake_dens = 0.
913	cv_gen = 0.
914	wake_cstar = 0.
915	ale_bl = 0.
916	ale_bl_trig = 0.
917	alp_bl = 0.
918	IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
919	IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
920	entr_therm = 0.
921	detr_therm = 0.
922	f0 = 0.
923	fm_therm = 0.
924	u_ancien(1, :) = u(:)
925	v_ancien(1, :) = v(:)
926
927	u10m = 0.
928	v10m = 0.
929	ale_wake = 0.
930	ale_bl_stat = 0.
931
932	!------------------------------------------------------------------------
933	! Make file containing restart for the physics (startphy.nc)
934
935	! NB: List of the variables to be written by phyredem (via put_field):
936	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
937	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
938	! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
939	! radsol,solsw,solswfdiff,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
940	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
941	! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
942	! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
943	! wake_deltat,wake_deltaq,wake_s,awake_s,wake_dens,awake_dens,cv_gen,wake_cstar,
944	! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
945
946	! NB2: The content of the startphy.nc file depends on some flags defined in
947	! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
948	! to be set at some arbitratry convenient values.
949	!------------------------------------------------------------------------
950	!Al1 =============== restart option ==========================
951	iflag_physiq = 0
952	CALL getin('iflag_physiq', iflag_physiq)
953
954	IF (.NOT.restart) THEN
955	iflag_pbl = 5
956	CALL phyredem ("startphy.nc")
957	else
958	! (desallocations)
959	PRINT*, 'callin surf final'
960	CALL pbl_surface_final(fder, snsrf, qsurfsrf, tsoil)
961	PRINT*, 'after surf final'
962	CALL fonte_neige_final(run_off_lic_0)
963	endif
964
965	ok_writedem = .FALSE.
966	PRINT*, 'apres phyredem'
967
968	END IF ! ok_writedem
969
970	!------------------------------------------------------------------------
971	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
972	! --------------------------------------------------
973	! NB: List of the variables to be written in limit.nc
974	! (by writelim.F, SUBROUTINE of 1DUTILS.h):
975	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
976	! phy_fter,phy_foce,phy_flic,phy_fsic)
977	!------------------------------------------------------------------------
978	DO i = 1, year_len
979	phy_nat(i) = nat_surf
980	phy_alb(i) = albedo
981	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
982	phy_rug(i) = rugos
983	phy_fter(i) = pctsrf(1, is_ter)
984	phy_foce(i) = pctsrf(1, is_oce)
985	phy_fsic(i) = pctsrf(1, is_sic)
986	phy_flic(i) = pctsrf(1, is_lic)
987	END DO
988
989	! fabrication de limit.nc
990	CALL writelim (1, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
991	phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
992
993	CALL phys_state_var_end
994	!Al1
995	IF (restart) THEN
996	PRINT*, 'CALL to restart dyn 1d'
997	Call dyn1deta0("start1dyn.nc", plev, play, phi, phis, presnivs, &
998	u, v, temp, q, omega2)
999
1000	PRINT*, 'fnday,annee_ref,day_ref,day_ini', &
1001	fnday, annee_ref, day_ref, day_ini
1002	!** CALL ymds2ju(annee_ref,mois,day_ini,heure,day)
1003	day = day_ini
1004	day_end = day_ini + nday
1005	daytime = day + time_ini / 24. ! 1st day and initial time of the simulation
1006
1007	! Print out the actual date of the beginning of the simulation :
1008	CALL ju2ymds(daytime, an, mois, jour, heure)
1009	PRINT *, ' Time of beginning : y m d h', an, mois, jour, heure / 3600.
1010
1011	day = int(daytime)
1012	time = daytime - day
1013
1014	PRINT, '*** intialised fields from restart1dyn *****'
1015	PRINT*, 'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
1016	PRINT*, 'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
1017	PRINT*, temp(1), q(1, 1), u(1), v(1), plev(1), phis
1018	! raz for safety
1019	DO l = 1, llm
1020	dq_dyn(l, 1) = 0.
1021	enddo
1022	END IF
1023	!Al1 ================ end restart =================================
1024	IF (ecrit_slab_oc==1) THEN
1025	open(97, file = 'div_slab.dat', STATUS = 'UNKNOWN')
1026	elseif (ecrit_slab_oc==0) THEN
1027	open(97, file = 'div_slab.dat', STATUS = 'OLD')
1028	END IF
1029
1030	!---------------------------------------------------------------------
1031	! Initialize target profile for RHT nudging if needed
1032	!---------------------------------------------------------------------
1033	IF (nudge(inudge_RHT)) THEN
1034	CALL nudge_RHT_init(plev, play, temp, q(:, 1), t_targ, rh_targ)
1035	END IF
1036	IF (nudge(inudge_UV)) THEN
1037	CALL nudge_UV_init(plev, play, u, v, u_targ, v_targ)
1038	END IF
1039
1040	!=====================================================================
1041	IF (CPPKEY_OUTPUTPHYSSCM) THEN
1042	CALL iophys_ini(timestep)
1043	END IF
1044	! START OF THE TEMPORAL LOOP :
1045	!=====================================================================
1046
1047	it_end = nint(fnday * day_step)
1048	!test JLD it_end = 10
1049	DO while(it<=it_end)
1050
1051	IF (prt_level>=1) THEN
1052	PRINT*, 'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
1053	it, day, time, it_end, day_step
1054	PRINT*, 'PAS DE TEMPS ', timestep
1055	endif
1056	!Al1 demande de restartphy.nc
1057	IF (it==it_end) lastcall = .True.
1058
1059	!---------------------------------------------------------------------
1060	! Geopotential :
1061	!---------------------------------------------------------------------
1062
1063	phi(1) = RD * temp(1) * (plev(1) - play(1)) / (.5 * (plev(1) + play(1)))
1064	DO l = 1, llm - 1
1065	phi(l + 1) = phi(l) + RD * (temp(l) + temp(l + 1)) * &
1066	(play(l) - play(l + 1)) / (play(l) + play(l + 1))
1067	enddo
1068
1069	!---------------------------------------------------------------------
1070	! Interpolation of forcings in time and onto model levels
1071	!---------------------------------------------------------------------
1072
1073	INCLUDE "old_1D_interp_cases.h"
1074
1075	IF (forcing_GCM2SCM) THEN
1076	write (, ) 'forcing_GCM2SCM not yet implemented'
1077	stop 'in time loop'
1078	END IF ! forcing_GCM2SCM
1079
1080	!!!!---------------------------------------------------------------------
1081	!!!! Geopotential :
1082	!!!!---------------------------------------------------------------------
1083	!!!
1084	!!! phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
1085	!!! do l = 1, llm-1
1086	!!! phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
1087	!!! & (play(l)-play(l+1))/(play(l)+play(l+1))
1088	!!! enddo
1089
1090	!---------------------------------------------------------------------
1091	! Listing output for debug prt_level>=1
1092	!---------------------------------------------------------------------
1093	IF (prt_level>=1) THEN
1094	PRINT *, ' avant physiq : -------- day time ', day, time
1095	WRITE(, ) 'firstcall,lastcall,phis', &
1096	firstcall, lastcall, phis
1097	end if
1098	IF (prt_level>=5) THEN
1099	WRITE(*, '(a10,2a4,4a13)') 'BEFOR1 IT=', 'it', 'l', &
1100	'presniv', 'plev', 'play', 'phi'
1101	WRITE(*, '(a10,2i4,4f13.2)') ('BEFOR1 IT= ', it, l, &
1102	presnivs(l), plev(l), play(l), phi(l), l = 1, llm)
1103	WRITE(*, '(a11,2a4,a11,6a8)') 'BEFOR2', 'it', 'l', &
1104	'presniv', 'u', 'v', 'temp', 'q1', 'q2', 'omega2'
1105	WRITE(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ', it, l, &
1106	presnivs(l), u(l), v(l), temp(l), q(l, 1), q(l, 2), omega2(l), l = 1, llm)
1107	END IF
1108
1109	!---------------------------------------------------------------------
1110	! Call physiq :
1111	!---------------------------------------------------------------------
1112	CALL physiq(ngrid, llm, &
1113	firstcall, lastcall, timestep, &
1114	plev, play, phi, phis, presnivs, &
1115	u, v, rot, temp, q, omega2, &
1116	du_phys, dv_phys, dt_phys, dq, dpsrf)
1117	firstcall = .FALSE.
1118
1119	!---------------------------------------------------------------------
1120	! Listing output for debug
1121	!---------------------------------------------------------------------
1122	IF (prt_level>=5) THEN
1123	WRITE(*, '(a11,2a4,4a13)') 'AFTER1 IT=', 'it', 'l', &
1124	'presniv', 'plev', 'play', 'phi'
1125	WRITE(*, '(a11,2i4,4f13.2)') ('AFTER1 it= ', it, l, &
1126	presnivs(l), plev(l), play(l), phi(l), l = 1, llm)
1127	WRITE(*, '(a11,2a4,a11,6a8)') 'AFTER2', 'it', 'l', &
1128	'presniv', 'u', 'v', 'temp', 'q1', 'q2', 'omega2'
1129	WRITE(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ', it, l, &
1130	presnivs(l), u(l), v(l), temp(l), q(l, 1), q(l, 2), omega2(l), l = 1, llm)
1131	WRITE(*, '(a11,2a4,a11,5a8)') 'AFTER3', 'it', 'l', &
1132	'presniv', 'du_phys', 'dv_phys', 'dt_phys', 'dq1', 'dq2'
1133	WRITE(*, '(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ', it, l, &
1134	presnivs(l), 86400 * du_phys(l), 86400 * dv_phys(l), &
1135	86400 * dt_phys(l), 86400 * dq(l, 1), dq(l, 2), l = 1, llm)
1136	WRITE(, ) 'dpsrf', dpsrf
1137	END IF
1138	!---------------------------------------------------------------------
1139	! Add physical tendencies :
1140	!---------------------------------------------------------------------
1141
1142	fcoriolis = 2. * sin(rpi * xlat / 180.) * romega
1143	IF (forcing_radconv .OR. forcing_fire) THEN
1144	fcoriolis = 0.0
1145	dt_cooling = 0.0
1146	d_t_adv = 0.0
1147	d_q_adv = 0.0
1148	END IF
1149	! PRINT*, 'calcul de fcoriolis ', fcoriolis
1150
1151	IF (forcing_toga .OR. forcing_GCSSold .OR. forcing_twpice &
1152	.OR.forcing_amma .OR. forcing_type==101) THEN
1153	fcoriolis = 0.0 ; ug = 0. ; vg = 0.
1154	END IF
1155
1156	IF(forcing_rico) THEN
1157	dt_cooling = 0.
1158	END IF
1159
1160	!CRio:Attention modif sp??cifique cas de Caroline
1161	IF (forcing_type==-1) THEN
1162	fcoriolis = 0.
1163	!Nudging
1164
1165	!on calcule dt_cooling
1166	DO l = 1, llm
1167	IF (play(l)>=20000.) THEN
1168	dt_cooling(l) = -1.5 / 86400.
1169	elseif ((play(l)>=10000.).AND.((play(l)<20000.))) THEN
1170	dt_cooling(l) = -1.5 / 86400. * (play(l) - 10000.) / (10000.) - 1. / 86400. * (20000. - play(l)) / 10000. * (temp(l) - 200.)
1171	else
1172	dt_cooling(l) = -1. * (temp(l) - 200.) / 86400.
1173	endif
1174	enddo
1175
1176	END IF
1177	!RC
1178	IF (forcing_sandu) THEN
1179	ug(1:llm) = u_mod(1:llm)
1180	vg(1:llm) = v_mod(1:llm)
1181	END IF
1182
1183	IF (prt_level >= 5) PRINT*, 'fcoriolis, xlat,mxcalc ', &
1184	fcoriolis, xlat, mxcalc
1185
1186	! PRINT *,'u-ug=',u-ug
1187
1188	!!!!!!!!!!!!!!!!!!!!!!!!
1189	! Geostrophic wind
1190	! Le calcul ci dessous est insuffisamment precis
1191	! du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
1192	! dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
1193	!!!!!!!!!!!!!!!!!!!!!!!!
1194	sfdt = sin(0.5 * fcoriolis * timestep)
1195	cfdt = cos(0.5 * fcoriolis * timestep)
1196	! PRINT *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
1197
1198	du_age(1:mxcalc) = -2. * sfdt / timestep * &
1199	(sfdt * (u(1:mxcalc) - ug(1:mxcalc)) - &
1200	cfdt * (v(1:mxcalc) - vg(1:mxcalc)))
1201	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
1202
1203	dv_age(1:mxcalc) = -2. * sfdt / timestep * &
1204	(cfdt * (u(1:mxcalc) - ug(1:mxcalc)) + &
1205	sfdt * (v(1:mxcalc) - vg(1:mxcalc)))
1206	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
1207
1208	!!!!!!!!!!!!!!!!!!!!!!!!
1209	! Nudging
1210	!!!!!!!!!!!!!!!!!!!!!!!!
1211	d_t_nudge(:) = 0.
1212	d_q_nudge(:, :) = 0.
1213	d_u_nudge(:) = 0.
1214	d_v_nudge(:) = 0.
1215	IF (nudge(inudge_RHT)) THEN
1216	CALL nudge_RHT(timestep, plev, play, t_targ, rh_targ, temp, q(:, 1), &
1217	d_t_nudge, d_q_nudge(:, 1))
1218	END IF
1219	IF (nudge(inudge_UV)) THEN
1220	CALL nudge_UV(timestep, plev, play, u_targ, v_targ, u, v, &
1221	d_u_nudge, d_v_nudge)
1222	END IF
1223
1224	IF (forcing_fire) THEN
1225
1226	!let ww=if ( alt le 1100 ) then alt*-0.00001 else 0
1227	!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt) else 0
1228	!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt) else 0
1229	d_t_adv = 0.
1230	d_q_adv = 0.
1231	teta = temp * (pzero / play)**rkappa
1232	d_t_adv = 0.
1233	d_q_adv = 0.
1234	DO l = 2, llm - 1
1235	IF (zlay(l)<=1100) THEN
1236	wwww = -0.00001 * zlay(l)
1237	d_t_adv(l) = -wwww * (teta(l) - teta(l + 1)) / (zlay(l) - zlay(l + 1)) / (pzero / play(l))**rkappa
1238	d_q_adv(l, 1:2) = -wwww * (q(l, 1:2) - q(l + 1, 1:2)) / (zlay(l) - zlay(l + 1))
1239	d_t_adv(l) = d_t_adv(l) + min(-3.75e-5, -7.5e-8 * zlay(l))
1240	d_q_adv(l, 1) = d_q_adv(l, 1) + max(1.5e-8, 3e-11 * zlay(l))
1241	endif
1242	enddo
1243
1244	END IF
1245
1246	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1247	! call writefield_phy('dv_age' ,dv_age,llm)
1248	! call writefield_phy('du_age' ,du_age,llm)
1249	! call writefield_phy('du_phys' ,du_phys,llm)
1250	! call writefield_phy('u_tend' ,u,llm)
1251	! call writefield_phy('u_g' ,ug,llm)
1252
1253	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1254	!! Increment state variables
1255	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1256	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
1257
1258	! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
1259	! au dessus de 700hpa, on relaxe vers les profils initiaux
1260	IF (forcing_sandu .OR. forcing_astex) THEN
1261	INCLUDE "1D_nudge_sandu_astex.h"
1262	else
1263	u(1:mxcalc) = u(1:mxcalc) + timestep * (&
1264	du_phys(1:mxcalc) &
1265	+ du_age(1:mxcalc) + du_adv(1:mxcalc) &
1266	+ d_u_nudge(1:mxcalc))
1267	v(1:mxcalc) = v(1:mxcalc) + timestep * (&
1268	dv_phys(1:mxcalc) &
1269	+ dv_age(1:mxcalc) + dv_adv(1:mxcalc) &
1270	+ d_v_nudge(1:mxcalc))
1271	q(1:mxcalc, :) = q(1:mxcalc, :) + timestep * (&
1272	dq(1:mxcalc, :) &
1273	+ d_q_adv(1:mxcalc, :) &
1274	+ d_q_nudge(1:mxcalc, :))
1275
1276	IF (prt_level>=3) THEN
1277	PRINT *, &
1278	'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ', &
1279	temp(1), dt_phys(1), d_t_adv(1), dt_cooling(1)
1280	PRINT*, 'dv_phys=', dv_phys
1281	PRINT*, 'dv_age=', dv_age
1282	PRINT*, 'dv_adv=', dv_adv
1283	PRINT*, 'd_v_nudge=', d_v_nudge
1284	PRINT*, v
1285	PRINT*, vg
1286	endif
1287
1288	temp(1:mxcalc) = temp(1:mxcalc) + timestep * (&
1289	dt_phys(1:mxcalc) &
1290	+ d_t_adv(1:mxcalc) &
1291	+ d_t_nudge(1:mxcalc) &
1292	+ dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
1293
1294	IF (CPPKEY_OUTPUTPHYSSCM) THEN
1295	CALL iophys_ecrit('d_t_adv', klev, 'd_t_adv', 'm/s', d_t_adv)
1296	CALL iophys_ecrit('d_t_nudge', klev, 'd_t_nudge', 'm/s', d_t_nudge)
1297	END IF
1298
1299	endif ! forcing_sandu or forcing_astex
1300
1301	teta = temp * (pzero / play)**rkappa
1302
1303	!---------------------------------------------------------------------
1304	! Nudge soil temperature if requested
1305	!---------------------------------------------------------------------
1306
1307	IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1308	ftsoil(1, isoil_nudge, :) = ftsoil(1, isoil_nudge, :) &
1309	- timestep / tau_soil_nudge * (ftsoil(1, isoil_nudge, :) - Tsoil_nudge)
1310	ENDIF
1311
1312	!---------------------------------------------------------------------
1313	! Add large-scale tendencies (advection, etc) :
1314	!---------------------------------------------------------------------
1315
1316	!cc nrlmd
1317	!cc tmpvar=teta
1318	!cc CALL advect_vert(llm,omega,timestep,tmpvar,plev)
1319	!cc
1320	!cc teta(1:mxcalc)=tmpvar(1:mxcalc)
1321	!cc tmpvar(:)=q(:,1)
1322	!cc CALL advect_vert(llm,omega,timestep,tmpvar,plev)
1323	!cc q(1:mxcalc,1)=tmpvar(1:mxcalc)
1324	!cc tmpvar(:)=q(:,2)
1325	!cc CALL advect_vert(llm,omega,timestep,tmpvar,plev)
1326	!cc q(1:mxcalc,2)=tmpvar(1:mxcalc)
1327
1328	END IF ! end if tendency of tendency should be added
1329
1330	!---------------------------------------------------------------------
1331	! Air temperature :
1332	!---------------------------------------------------------------------
1333	IF (lastcall) THEN
1334	PRINT*, 'Pas de temps final ', it
1335	CALL ju2ymds(daytime, an, mois, jour, heure)
1336	PRINT*, 'a la date : a m j h', an, mois, jour, heure / 3600.
1337	END IF
1338
1339	! incremente day time
1340	! PRINT*,'daytime bef',daytime,1./day_step
1341	daytime = daytime + 1. / day_step
1342	!Al1dbg
1343	day = int(daytime + 0.1 / day_step)
1344	! time = max(daytime-day,0.0)
1345	!Al1&jyg: correction de bug
1346	!cc time = real(mod(it,day_step))/day_step
1347	time = time_ini / 24. + real(mod(it, day_step)) / day_step
1348	! PRINT*,'daytime nxt time',daytime,time
1349	it = it + 1
1350
1351	END DO
1352
1353	!Al1
1354	IF (ecrit_slab_oc/=-1) close(97)
1355
1356	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1357	! -------------------------------------
1358	CALL dyn1dredem("restart1dyn.nc", &
1359	plev, play, phi, phis, presnivs, &
1360	u, v, temp, q, omega2)
1361
1362	CALL abort_gcm ('lmdz1d ', 'The End ', 0)
1363
1364	END SUBROUTINE old_lmdz1d
1365
1366	INCLUDE "old_1DUTILS_read_interp.h"
1367	END MODULE lmdz_old_lmdz1d

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

Download in other formats:

Original Format