Context Navigation

lmdz_old_lmdz1d.F90 @ 5104

Last change on this file since 5104 was 5104, checked in by abarral, 2 months ago
Replace 1DUTILS.h by module lmdz_1dutils.f90 Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files) Convert .F in DUST to .f90
File size: 54.9 KB

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

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

Download in other formats:

Original Format