Context Navigation

old_lmdz1d.f90 @ 5601

Last change on this file since 5601 was 5601, checked in by aborella, 2 months ago

Multiple changes:

tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

added a new tracer for total specific humidity in contrails
rework of the mixing process for cirrus clouds (and contrails)
changed the numerical integration of ice crystals' sublimation
subroutines do not take real inputs anymore (at least klon tables)
added more radiative diagnostics for contrails

File size: 51.5 KB

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

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

Download in other formats:

Original Format