Context Navigation

lmdz1d.F @ 1943

Last change on this file since 1943 was 1943, checked in by idelkadi, 11 years ago

Changes concerinng the Hourdin et al. 2002 version of the thermals and
Ayotte cases

calltherm.F90

call to thermcell_2002 modified :
1) iflag_thermals changed from 1 to >= 1000 ; iflag_thermals-1000
controls sub options.
2) thermals w and fraction in output files.

hbtm.F

Singularity in the 1/heat dependency of the Monin Obukov length L when
heat=0. Since 1/L is used rather than L, it is preferable to compute
directly L. There is a dependency in 1/u* then which is treated with a
threshold value.
(+ some cleaning in the syntax).

lmdz1d.F

If nday<0, -nday is used as the total number of time steps of the
simulation.
The option with imposed wtsurf and wqsurf read in lmdz1d.def was not
active anymore.
< IF(.NOT.ok_flux_surf) THEN
changed to

IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN

before

fsens=-wtsurf*rcpd*rho(1)
flat=-wqsurf*rlvtt*rho(1)

phys_output_write.F90 et phys_local_var_mod.F90

Removing the d_u_ajsb contribution to duthe (same for dv).
Those tendencies are not computed by the model ...
< zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
< d_u_ajsb(1:klon,1:klev)/pdtphys
---

zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys

! d_u_ajsb(1:klon,1:klev)/pdtphys

thermcell_dq.F90, thermcell_main.F90

Some cleaning

thermcell_old.F

Control by iflag_thermals >= 1000
wa_moy and fraca in outputs
+ cleaning

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 35.0 KB

Line
1	PROGRAM lmdz1d
2
3	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
4	use phys_state_var_mod
5	use comgeomphy
6	use dimphy
7	use surface_data, only : type_ocean,ok_veget
8	use pbl_surface_mod, only : ftsoil, pbl_surface_init,
9	$ pbl_surface_final
10	use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
11
12	use infotrac ! new
13	use control_mod
14	USE indice_sol_mod
15	USE phyaqua_mod
16
17	implicit none
18	#include "dimensions.h"
19	#include "YOMCST.h"
20	#include "temps.h"
21	!!#include "control.h"
22	#include "iniprint.h"
23	#include "clesphys.h"
24	#include "dimsoil.h"
25	!#include "indicesol.h"
26
27	#include "comvert.h"
28	#include "compar1d.h"
29	#include "flux_arp.h"
30	#include "tsoilnudge.h"
31	#include "fcg_gcssold.h"
32	!!!#include "fbforcing.h"
33
34	!=====================================================================
35	! DECLARATIONS
36	!=====================================================================
37
38	!---------------------------------------------------------------------
39	! Externals
40	!---------------------------------------------------------------------
41	external fq_sat
42	real fq_sat
43
44	!---------------------------------------------------------------------
45	! Arguments d' initialisations de la physique (USER DEFINE)
46	!---------------------------------------------------------------------
47
48	integer, parameter :: ngrid=1
49	real :: zcufi = 1.
50	real :: zcvfi = 1.
51
52	!- real :: nat_surf
53	!- logical :: ok_flux_surf
54	!- real :: fsens
55	!- real :: flat
56	!- real :: tsurf
57	!- real :: rugos
58	!- real :: qsol(1:2)
59	!- real :: qsurf
60	!- real :: psurf
61	!- real :: zsurf
62	!- real :: albedo
63	!-
64	!- real :: time = 0.
65	!- real :: time_ini
66	!- real :: xlat
67	!- real :: xlon
68	!- real :: wtsurf
69	!- real :: wqsurf
70	!- real :: restart_runoff
71	!- real :: xagesno
72	!- real :: qsolinp
73	!- real :: zpicinp
74	!-
75	real :: fnday
76	real :: day, daytime
77	real :: day1
78	real :: heure
79	integer :: jour
80	integer :: mois
81	integer :: an
82
83	!
84	real :: paire = 1. ! aire de la maille
85	!** common /flux_arp/fsens,flat,ok_flux_surf
86
87	!---------------------------------------------------------------------
88	! Declarations related to forcing and initial profiles
89	!---------------------------------------------------------------------
90
91	integer :: kmax = llm
92	integer nlev_max,llm700
93	parameter (nlev_max = 1000)
94	real timestep, frac, timeit
95	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max),
96	. uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
97	. ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
98	. dqtdxls(nlev_max),dqtdyls(nlev_max),
99	. dqtdtls(nlev_max),thlpcar(nlev_max)
100
101	real :: fff
102	c integer :: forcing_type
103	logical :: forcing_les = .false.
104	logical :: forcing_armcu = .false.
105	logical :: forcing_rico = .false.
106	logical :: forcing_radconv = .false.
107	logical :: forcing_toga = .false.
108	logical :: forcing_twpice = .false.
109	logical :: forcing_amma = .false.
110	logical :: forcing_GCM2SCM = .false.
111	logical :: forcing_GCSSold = .false.
112	logical :: forcing_sandu = .false.
113	logical :: forcing_astex = .false.
114	logical :: forcing_fire = .false.
115	integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
116	! (cf read_tsurf1d.F)
117
118	!vertical advection computation
119	! real d_t_z(llm), d_q_z(llm)
120	! real d_t_dyn_z(llm), d_q_dyn_z(llm)
121	! real zz(llm)
122	! real zfact
123
124	!flag forcings
125	logical :: nudge_wind=.true.
126	logical :: nudge_thermo=.false.
127	logical :: cptadvw=.true.
128	!=====================================================================
129	! DECLARATIONS FOR EACH CASE
130	!=====================================================================
131	!
132	#include "1D_decl_cases.h"
133	!
134	!---------------------------------------------------------------------
135	! Declarations related to vertical discretization:
136	!---------------------------------------------------------------------
137	real :: pzero=1.e5
138	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
139	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
140
141	!---------------------------------------------------------------------
142	! Declarations related to variables
143	!---------------------------------------------------------------------
144
145	integer :: iq
146	real :: phi(llm)
147	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
148	real :: rlat_rad(1),rlon_rad(1)
149	real :: omega(llm+1),omega2(llm),rho(llm+1)
150	real :: ug(llm),vg(llm),fcoriolis
151	real :: sfdt, cfdt
152	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
153	real :: du_dyn(llm),dv_dyn(llm),dt_dyn(llm)
154	real :: dt_cooling(llm),d_t_cool(llm),d_th_adv(llm)
155	real :: dq_cooling(llm),d_q_cool(llm)
156	real :: tmpvar(llm)
157	real :: alpha
158
159	REAL, ALLOCATABLE, DIMENSION(:,:):: q
160	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
161	REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
162	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
163
164	!---------------------------------------------------------------------
165	! Initialization of surface variables
166	!---------------------------------------------------------------------
167	real :: run_off_lic_0(1)
168	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
169	real :: evap(1,nbsrf),frugs(1,nbsrf)
170	real :: tsoil(1,nsoilmx,nbsrf)
171	real :: agesno(1,nbsrf)
172
173	!---------------------------------------------------------------------
174	! Call to phyredem
175	!---------------------------------------------------------------------
176	logical :: ok_writedem =.true.
177
178	!---------------------------------------------------------------------
179	! Call to physiq
180	!---------------------------------------------------------------------
181	integer, parameter :: longcles=20
182	logical :: firstcall=.true.
183	logical :: lastcall=.false.
184	real :: phis = 0.0
185	real :: clesphy0(longcles) = 0.0
186	real :: dpsrf
187
188	!---------------------------------------------------------------------
189	! Initializations of boundary conditions
190	!---------------------------------------------------------------------
191	integer, parameter :: yd = 360
192	real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
193	real :: phy_alb (yd) ! Albedo land only (old value condsurf_jyg=0.3)
194	real :: phy_sst (yd) ! SST (will not be used; cf read_tsurf1d.F)
195	real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
196	real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
197	real :: phy_ice (yd) = 0.0 ! Fraction de glace
198	real :: phy_fter(yd) = 0.0 ! Fraction de terre
199	real :: phy_foce(yd) = 0.0 ! Fraction de ocean
200	real :: phy_fsic(yd) = 0.0 ! Fraction de glace
201	real :: phy_flic(yd) = 0.0 ! Fraction de glace
202
203	!---------------------------------------------------------------------
204	! Fichiers et d'autres variables
205	!---------------------------------------------------------------------
206	real ttt,bow,q1
207	integer :: ierr,k,l,i,it=1,mxcalc
208	integer jjmp1
209	parameter (jjmp1=jjm+1-1/jjm)
210	INTEGER nbteta
211	PARAMETER(nbteta=3)
212	REAL dudyn(iim+1,jjmp1,llm)
213	REAL PVteta(1,nbteta)
214	INTEGER read_climoz
215	!Al1
216	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
217	data ecrit_slab_oc/-1/
218
219	!=====================================================================
220	! INITIALIZATIONS
221	!=====================================================================
222	! Initialization of Common turb_forcing
223	dtime_frcg = 0.
224	Turb_fcg_gcssold=.false.
225	hthturb_gcssold = 0.
226	hqturb_gcssold = 0.
227
228	!---------------------------------------------------------------------
229	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
230	!---------------------------------------------------------------------
231	cAl1
232	call conf_unicol(99)
233	cAl1 moves this gcssold var from common fcg_gcssold to
234	Turb_fcg_gcssold = xTurb_fcg_gcssold
235	c --------------------------------------------------------------------
236	close(1)
237	cAl1
238	write(,) 'lmdz1d.def lu => unicol.def'
239
240	! forcing_type defines the way the SCM is forced:
241	!forcing_type = 0 ==> forcing_les = .true.
242	! initial profiles from file prof.inp.001
243	! no forcing by LS convergence ;
244	! surface temperature imposed ;
245	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
246	!forcing_type = 1 ==> forcing_radconv = .true.
247	! idem forcing_type = 0, but the imposed radiative cooling
248	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
249	! then there is no radiative cooling at all)
250	!forcing_type = 2 ==> forcing_toga = .true.
251	! initial profiles from TOGA-COARE IFA files
252	! LS convergence and SST imposed from TOGA-COARE IFA files
253	!forcing_type = 3 ==> forcing_GCM2SCM = .true.
254	! initial profiles from the GCM output
255	! LS convergence imposed from the GCM output
256	!forcing_type = 4 ==> forcing_twpice = .true.
257	! initial profiles from TWP-ICE cdf file
258	! LS convergence, omega and SST imposed from TWP-ICE files
259	!forcing_type = 5 ==> forcing_rico = .true.
260	! initial profiles from RICO files
261	! LS convergence imposed from RICO files
262	!forcing_type = 6 ==> forcing_amma = .true.
263	! initial profiles from AMMA nc file
264	! LS convergence, omega and surface fluxes imposed from AMMA file
265	!forcing_type = 40 ==> forcing_GCSSold = .true.
266	! initial profile from GCSS file
267	! LS convergence imposed from GCSS file
268	!forcing_type = 50 ==> forcing_fire = .true.
269	! forcing from fire.nc
270	!forcing_type = 59 ==> forcing_sandu = .true.
271	! initial profiles from sanduref file: see prof.inp.001
272	! SST varying with time and divergence constante: see ifa_sanduref.txt file
273	! Radiation has to be computed interactively
274	!forcing_type = 60 ==> forcing_astex = .true.
275	! initial profiles from file: see prof.inp.001
276	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
277	! Radiation has to be computed interactively
278	!forcing_type = 61 ==> forcing_armcu = .true.
279	! initial profiles from file: see prof.inp.001
280	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
281	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
282	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
283	! Radiation to be switched off
284	!
285	if (forcing_type .eq.0) THEN
286	forcing_les = .true.
287	elseif (forcing_type .eq.1) THEN
288	forcing_radconv = .true.
289	elseif (forcing_type .eq.2) THEN
290	forcing_toga = .true.
291	elseif (forcing_type .eq.3) THEN
292	forcing_GCM2SCM = .true.
293	elseif (forcing_type .eq.4) THEN
294	forcing_twpice = .true.
295	elseif (forcing_type .eq.5) THEN
296	forcing_rico = .true.
297	elseif (forcing_type .eq.6) THEN
298	forcing_amma = .true.
299	elseif (forcing_type .eq.40) THEN
300	forcing_GCSSold = .true.
301	elseif (forcing_type .eq.50) THEN
302	forcing_fire = .true.
303	elseif (forcing_type .eq.59) THEN
304	forcing_sandu = .true.
305	elseif (forcing_type .eq.60) THEN
306	forcing_astex = .true.
307	elseif (forcing_type .eq.61) THEN
308	forcing_armcu = .true.
309	IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
310	else
311	write (,) 'ERROR : unknown forcing_type ', forcing_type
312	stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
313	ENDIF
314	print*,"forcing type=",forcing_type
315
316	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
317	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
318	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
319	! through the common sst_forcing.
320
321	type_ts_forcing = 0
322	if (forcing_toga .or. forcing_sandu .or. forcing_astex)
323	: type_ts_forcing = 1
324
325	!---------------------------------------------------------------------
326	! Definition of the run
327	!---------------------------------------------------------------------
328
329	call conf_gcm( 99, .TRUE. , clesphy0 )
330	c-----------------------------------------------------------------------
331	c Choix du calendrier
332	c -------------------
333
334	c calend = 'earth_365d'
335	if (calend == 'earth_360d') then
336	call ioconf_calendar('360d')
337	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
338	else if (calend == 'earth_365d') then
339	call ioconf_calendar('noleap')
340	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
341	else if (calend == 'earth_366d') then
342	call ioconf_calendar('all_leap')
343	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
344	else if (calend == 'gregorian') then
345	call ioconf_calendar('gregorian') ! not to be used by normal users
346	write(,)'CALENDRIER CHOISI: Gregorien'
347	else
348	write (,) 'ERROR : unknown calendar ', calend
349	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
350	endif
351	c-----------------------------------------------------------------------
352	c
353	cc Date :
354	c La date est supposee donnee sous la forme [annee, numero du jour dans
355	c l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
356	c On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
357	c Le numero du jour est dans "day". L heure est traitee separement.
358	c La date complete est dans "daytime" (l'unite est le jour).
359	if (nday>0) then
360	fnday=nday
361	else
362	fnday=-nday/float(day_step)
363	endif
364
365	c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
366	IF(forcing_type .EQ. 61) fnday=53100./86400.
367	c Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
368	IF(forcing_type .EQ. 6) fnday=64800./86400.
369	annee_ref = anneeref
370	mois = 1
371	day_ref = dayref
372	heure = 0.
373	itau_dyn = 0
374	itau_phy = 0
375	call ymds2ju(annee_ref,mois,day_ref,heure,day)
376	day_ini = day
377	day_end = day_ini + fnday
378
379	IF (forcing_type .eq.2) THEN
380	! Convert the initial date of Toga-Coare to Julian day
381	call ymds2ju
382	$ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
383
384	ELSEIF (forcing_type .eq.4) THEN
385	! Convert the initial date of TWPICE to Julian day
386	call ymds2ju
387	$ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
388	$ ,day_ju_ini_twpi)
389	ELSEIF (forcing_type .eq.6) THEN
390	! Convert the initial date of AMMA to Julian day
391	call ymds2ju
392	$ (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma
393	$ ,day_ju_ini_amma)
394
395	ELSEIF (forcing_type .eq.59) THEN
396	! Convert the initial date of Sandu case to Julian day
397	call ymds2ju
398	$ (year_ini_sandu,mth_ini_sandu,day_ini_sandu,
399	$ time_ini*3600.,day_ju_ini_sandu)
400
401	ELSEIF (forcing_type .eq.60) THEN
402	! Convert the initial date of Astex case to Julian day
403	call ymds2ju
404	$ (year_ini_astex,mth_ini_astex,day_ini_astex,
405	$ time_ini*3600.,day_ju_ini_astex)
406
407	ELSEIF (forcing_type .eq.61) THEN
408
409	! Convert the initial date of Arm_cu case to Julian day
410	call ymds2ju
411	$ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
412	$ ,day_ju_ini_armcu)
413	ENDIF
414
415	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
416	! Print out the actual date of the beginning of the simulation :
417	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
418	print *,' Time of beginning : ',
419	$ year_print, month_print, day_print, sec_print
420
421	!---------------------------------------------------------------------
422	! Initialization of dimensions, geometry and initial state
423	!---------------------------------------------------------------------
424	call init_phys_lmdz(1,1,llm,1,(/1/))
425	call suphel
426	call initcomgeomphy
427	call infotrac_init
428
429	allocate(q(llm,nqtot)) ; q(:,:)=0.
430	allocate(dq(llm,nqtot))
431	allocate(dq_dyn(llm,nqtot))
432	allocate(d_q_adv(llm,nqtot))
433
434	c
435	c No ozone climatology need be read in this pre-initialization
436	c (phys_state_var_init is called again in physiq)
437	read_climoz = 0
438	c
439	call phys_state_var_init(read_climoz)
440
441	if (ngrid.ne.klon) then
442	print*,'stop in inifis'
443	print*,'Probleme de dimensions :'
444	print*,'ngrid = ',ngrid
445	print*,'klon = ',klon
446	stop
447	endif
448	!!!=====================================================================
449	!!! Feedback forcing values for Gateaux differentiation (al1)
450	!!!=====================================================================
451	!!! Surface Planck forcing bracketing call radiation
452	!! surf_Planck = 0.
453	!! surf_Conv = 0.
454	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
455	!!! a mettre dans le lmdz1d.def ou autre
456	!!
457	!!
458	qsol = qsolinp
459	qsurf = fq_sat(tsurf,psurf/100.)
460	rlat=xlat
461	rlon=xlon
462	day1= day_ini
463	time=daytime-day
464	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
465	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
466
467	!
468	!! mpl et jyg le 22/08/2012 :
469	!! pour que les cas a flux de surface imposes marchent
470	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
471	fsens=-wtsurfrcpdrho(1)
472	flat=-wqsurfrlvttrho(1)
473	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
474	ENDIF
475	print*,'Flux sol ',fsens,flat
476	!! ok_flux_surf=.false.
477	!! fsens=-wtsurfrcpdrho(1)
478	!! flat=-wqsurfrlvttrho(1)
479	!!!!
480
481	! Vertical discretization and pressure levels at half and mid levels:
482
483	pa = 5e4
484	!! preff= 1.01325e5
485	preff = psurf
486	IF (ok_old_disvert) THEN
487	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
488	print *,'On utilise disvert0'
489	ELSE
490	call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
491	: scaleheight)
492	print *,'On utilise disvert'
493	c Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
494	c Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
495	ENDIF
496	sig_s=presnivs/preff
497	plev =ap+bp*psurf
498	play = 0.5*(plev(1:llm)+plev(2:llm+1))
499	ccc zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles
500
501	IF (forcing_type .eq. 59) THEN
502	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
503	write(,) '***********************'
504	do l = 1, llm
505	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
506	if (trouve_700 .and. play(l).le.70000) then
507	llm700=l
508	print *,'llm700,play=',llm700,play(l)/100.
509	trouve_700= .false.
510	endif
511	enddo
512	write(,) '***********************'
513	ENDIF
514
515	c
516	!=====================================================================
517	! EVENTUALLY, READ FORCING DATA :
518	!=====================================================================
519
520	#include "1D_read_forc_cases.h"
521
522	if (forcing_GCM2SCM) then
523	write (,) 'forcing_GCM2SCM not yet implemented'
524	stop 'in initialization'
525	endif ! forcing_GCM2SCM
526
527	print*,'mxcalc=',mxcalc
528	print*,'zlay=',zlay(mxcalc)
529	print*,'play=',play(mxcalc)
530
531	cAl1 pour SST forced, appellé depuis ocean_forced_noice
532	ts_cur = tsurf ! SST used in read_tsurf1d
533	!=====================================================================
534	! Initialisation de la physique :
535	!=====================================================================
536
537	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
538	!
539	! day_step, iphysiq lus dans gcm.def ci-dessus
540	! timestep: calcule ci-dessous from rday et day_step
541	! ngrid=1
542	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
543	! rday: defini dans suphel.F (86400.)
544	! day_ini: lu dans run.def (dayref)
545	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
546	! airefi,zcufi,zcvfi initialises au debut de ce programme
547	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
548	day_step = float(nsplit_phys)*day_step/float(iphysiq)
549	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
550	timestep =rday/day_step
551	dtime_frcg = timestep
552	!
553	zcufi=airefi
554	zcvfi=airefi
555	!
556	rlat_rad(:)=rlat(:)*rpi/180.
557	rlon_rad(:)=rlon(:)*rpi/180.
558
559	call iniphysiq(ngrid,llm,rday,day_ini,timestep,
560	. rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
561	print*,'apres iniphysiq'
562
563	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
564	co2_ppm= 330.0
565	solaire=1370.0
566
567	! Ecriture du startphy avant le premier appel a la physique.
568	! On le met juste avant pour avoir acces a tous les champs
569	! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
570
571	if (ok_writedem) then
572
573	!--------------------------------------------------------------------------
574	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
575	! need : qsol fder snow qsurf evap rugos agesno ftsoil
576	!--------------------------------------------------------------------------
577
578	type_ocean = "force"
579	run_off_lic_0(1) = restart_runoff
580	call fonte_neige_init(run_off_lic_0)
581
582	fder=0.
583	snsrf(1,:)=0. ! couverture de neige des sous surface
584	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
585	evap=0.
586	frugs(1,:)=rugos ! couverture de neige des sous surface
587	agesno = xagesno
588	tsoil(:,:,:)=tsurf
589	!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
590	! tsoil(1,1,1)=299.18
591	! tsoil(1,2,1)=300.08
592	! tsoil(1,3,1)=301.88
593	! tsoil(1,4,1)=305.48
594	! tsoil(1,5,1)=308.00
595	! tsoil(1,6,1)=308.00
596	! tsoil(1,7,1)=308.00
597	! tsoil(1,8,1)=308.00
598	! tsoil(1,9,1)=308.00
599	! tsoil(1,10,1)=308.00
600	! tsoil(1,11,1)=308.00
601	!-----------------------------------------------------------------------
602	call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,
603	& evap, frugs, agesno, tsoil)
604
605	!------------------ prepare limit conditions for limit.nc -----------------
606	!-- Ocean force
607
608	print*,'avant phyredem'
609	pctsrf(1,:)=0.
610	if (nat_surf.eq.0.) then
611	pctsrf(1,is_oce)=1.
612	pctsrf(1,is_ter)=0.
613	else
614	pctsrf(1,is_oce)=0.
615	pctsrf(1,is_ter)=1.
616	end if
617
618	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf
619	$ ,pctsrf(1,is_oce),pctsrf(1,is_ter)
620
621	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
622	zpic = zpicinp
623	ftsol=tsurf
624	falb1 = albedo
625	falb2 = albedo
626	rugoro=rugos
627	t_ancien(1,:)=temp(:)
628	q_ancien(1,:)=q(:,1)
629	pbl_tke=1.e-8
630
631	rain_fall=0.
632	snow_fall=0.
633	solsw=0.
634	sollw=0.
635	radsol=0.
636	rnebcon=0.
637	ratqs=0.
638	clwcon=0.
639	zmea=0.
640	zstd=0.
641	zsig=0.
642	zgam=0.
643	zval=0.
644	zthe=0.
645	sig1=0.
646	w01=0.
647	u_ancien(1,:)=u(:)
648	v_ancien(1,:)=v(:)
649
650	!------------------------------------------------------------------------
651	! Make file containing restart for the physics (startphy.nc)
652	!
653	! NB: List of the variables to be written by phyredem (via put_field):
654	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
655	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
656	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
657	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
658	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
659	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
660	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,sig1,w01
661	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
662	!------------------------------------------------------------------------
663	CAl1 =============== restart option ==========================
664	if (.not.restart) then
665	call phyredem ("startphy.nc")
666	else
667	c (desallocations)
668	print*,'callin surf final'
669	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,
670	& evap, frugs, agesno, tsoil)
671	print*,'after surf final'
672	CALL fonte_neige_final(run_off_lic_0)
673	endif
674
675	ok_writedem=.false.
676	print*,'apres phyredem'
677
678	endif ! ok_writedem
679
680	!------------------------------------------------------------------------
681	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
682	! --------------------------------------------------
683	! NB: List of the variables to be written in limit.nc
684	! (by writelim.F, subroutine of 1DUTILS.h):
685	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
686	! phy_fter,phy_foce,phy_flic,phy_fsic)
687	!------------------------------------------------------------------------
688	do i=1,yd
689	phy_nat(i) = nat_surf
690	phy_alb(i) = albedo
691	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
692	phy_rug(i) = rugos
693	phy_fter(i) = pctsrf(1,is_ter)
694	phy_foce(i) = pctsrf(1,is_oce)
695	phy_fsic(i) = pctsrf(1,is_sic)
696	phy_flic(i) = pctsrf(1,is_lic)
697	enddo
698
699	C fabrication de limit.nc
700	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,
701	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
702
703
704	call phys_state_var_end
705	cAl1
706	if (restart) then
707	print*,'call to restart dyn 1d'
708	Call dyn1deta0("start1dyn.nc",
709	& plev,play,phi,phis,presnivs,
710	& u,v,temp,q,omega2)
711
712	print*,'fnday,annee_ref,day_ref,day_ini',
713	& fnday,annee_ref,day_ref,day_ini
714	c** call ymds2ju(annee_ref,mois,day_ini,heure,day)
715	day = day_ini
716	day_end = day_ini + nday
717	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
718
719	! Print out the actual date of the beginning of the simulation :
720	call ju2ymds(daytime, an, mois, jour, heure)
721	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
722
723	day = int(daytime)
724	time=daytime-day
725
726	print,'*** intialised fields from restart1dyn *****'
727	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
728	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
729	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
730	c raz for safety
731	do l=1,llm
732	dq_dyn(l,1) = 0.
733	enddo
734	endif
735	!Al1 ================ end restart =================================
736	IF (ecrit_slab_oc.eq.1) then
737	open(97,file='div_slab.dat',STATUS='UNKNOWN')
738	elseif (ecrit_slab_oc.eq.0) then
739	open(97,file='div_slab.dat',STATUS='OLD')
740	endif
741	!=====================================================================
742	! START OF THE TEMPORAL LOOP :
743	!=====================================================================
744
745	do while(it.le.nint(fnday*day_step))
746
747	if (prt_level.ge.1) then
748	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',
749	. it,day,time,nint(fnday*day_step)
750	print*,'PAS DE TEMPS ',timestep
751	endif
752	!Al1 demande de restartphy.nc
753	if (it.eq.nint(fnday*day_step)) lastcall=.True.
754
755	!---------------------------------------------------------------------
756	! Interpolation of forcings in time and onto model levels
757	!---------------------------------------------------------------------
758
759	#include "1D_interp_cases.h"
760
761	if (forcing_GCM2SCM) then
762	write (,) 'forcing_GCM2SCM not yet implemented'
763	stop 'in time loop'
764	endif ! forcing_GCM2SCM
765
766	!---------------------------------------------------------------------
767	! Geopotential :
768	!---------------------------------------------------------------------
769
770	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
771	do l = 1, llm-1
772	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1))
773	. (play(l)-play(l+1))/(play(l)+play(l+1))
774	enddo
775
776	!---------------------------------------------------------------------
777	! Listing output for debug prt_level>=1
778	!---------------------------------------------------------------------
779	if (prt_level>=1) then
780	print *,' avant physiq : -------- day time ',day,time
781	write(,) 'firstcall,lastcall,phis',
782	: firstcall,lastcall,phis
783	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',
784	: 'presniv','plev','play','phi'
785	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,
786	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
787	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',
788	: 'presniv','u','v','temp','q1','q2','omega2'
789	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,
790	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
791	endif
792
793	!---------------------------------------------------------------------
794	! Call physiq :
795	!---------------------------------------------------------------------
796
797	call physiq(ngrid,llm,
798	: firstcall,lastcall,
799	: day,time,timestep,
800	: plev,play,phi,phis,presnivs,clesphy0,
801	: u,v,temp,q,omega2,
802	: du_phys,dv_phys,dt_phys,dq,dpsrf,
803	: dudyn,PVteta)
804	firstcall=.false.
805
806	!---------------------------------------------------------------------
807	! Listing output for debug prt_level>=1
808	!---------------------------------------------------------------------
809	if (prt_level>=1) then
810	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',
811	: 'presniv','plev','play','phi'
812	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,
813	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
814	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',
815	: 'presniv','u','v','temp','q1','q2','omega2'
816	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,
817	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
818	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',
819	: 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
820	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,
821	: presnivs(l),86400du_phys(l),86400dv_phys(l),
822	: 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
823	write(,) 'dpsrf',dpsrf
824	endif
825	!---------------------------------------------------------------------
826	! Add physical tendencies :
827	!---------------------------------------------------------------------
828
829	fcoriolis=2.sin(rpixlat/180.)*romega
830	if (forcing_radconv .or. forcing_fire) then
831	fcoriolis=0.0
832	dt_cooling=0.0
833	d_th_adv=0.0
834	d_q_adv=0.0
835	endif
836
837	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice
838	: .or.forcing_amma) then
839	fcoriolis=0.0 ; ug=0. ; vg=0.
840	endif
841	if(forcing_rico) then
842	dt_cooling=0.
843	endif
844
845	print*, 'fcoriolis ', fcoriolis, xlat
846
847	du_age(1:mxcalc)=
848	: fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
849	dv_age(1:mxcalc)=
850	: -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
851
852	!!!!!!!!!!!!!!!!!!!!!!!!
853	! Geostrophic wind
854	!!!!!!!!!!!!!!!!!!!!!!!!
855	sfdt = sin(0.5fcoriolistimestep)
856	cfdt = cos(0.5fcoriolistimestep)
857	!
858	du_age(1:mxcalc)= -2.sfdt/timestep
859	: (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -
860	: cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
861	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
862	!
863	dv_age(1:mxcalc)= -2.sfdt/timestep
864	: (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +
865	: sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
866	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
867	!
868	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
869	! call writefield_phy('dv_age' ,dv_age,llm)
870	! call writefield_phy('du_age' ,du_age,llm)
871	! call writefield_phy('du_phys' ,du_phys,llm)
872	! call writefield_phy('u_tend' ,u,llm)
873	! call writefield_phy('u_g' ,ug,llm)
874	!
875	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
876	!! Increment state variables
877	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
878	u(1:mxcalc)=u(1:mxcalc) + timestep*(
879	: du_phys(1:mxcalc)
880	: +du_age(1:mxcalc) )
881	v(1:mxcalc)=v(1:mxcalc) + timestep*(
882	: dv_phys(1:mxcalc)
883	: +dv_age(1:mxcalc) )
884	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
885	: dq(1:mxcalc,:)
886	: +d_q_adv(1:mxcalc,:) )
887
888	if (prt_level.ge.1) then
889	print *,
890	: 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',
891	: temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
892	print*,du_phys
893	print*, v
894	print*, vg
895	endif
896
897	temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
898	. dt_phys(1:mxcalc)
899	. +d_th_adv(1:mxcalc)
900	. +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
901
902	teta=temp(pzero/play)*rkappa
903	!
904	!---------------------------------------------------------------------
905	! Nudge soil temperature if requested
906	!---------------------------------------------------------------------
907
908	IF (nudge_tsoil) THEN
909	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)
910	. -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
911	ENDIF
912
913	!---------------------------------------------------------------------
914	! Add large-scale tendencies (advection, etc) :
915	!---------------------------------------------------------------------
916
917	ccc nrlmd
918	ccc tmpvar=teta
919	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
920	ccc
921	ccc teta(1:mxcalc)=tmpvar(1:mxcalc)
922	ccc tmpvar(:)=q(:,1)
923	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
924	ccc q(1:mxcalc,1)=tmpvar(1:mxcalc)
925	ccc tmpvar(:)=q(:,2)
926	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
927	ccc q(1:mxcalc,2)=tmpvar(1:mxcalc)
928
929	!---------------------------------------------------------------------
930	! Air temperature :
931	!---------------------------------------------------------------------
932	if (lastcall) then
933	print*,'Pas de temps final ',it
934	call ju2ymds(daytime, an, mois, jour, heure)
935	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
936	endif
937
938	c incremente day time
939	c print*,'daytime bef',daytime,1./day_step
940	daytime = daytime+1./day_step
941	!Al1dbg
942	day = int(daytime+0.1/day_step)
943	c time = max(daytime-day,0.0)
944	cAl1&jyg: correction de bug
945	ccc time = real(mod(it,day_step))/day_step
946	time = time_ini/24.+real(mod(it,day_step))/day_step
947	c print*,'daytime nxt time',daytime,time
948	it=it+1
949
950	enddo
951
952	!Al1
953	if (ecrit_slab_oc.ne.-1) close(97)
954
955	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
956	! -------------------------------------
957	call dyn1dredem("restart1dyn.nc",
958	: plev,play,phi,phis,presnivs,
959	: u,v,temp,q,omega2)
960
961	CALL abort_gcm ('lmdz1d ','The End ',0)
962
963	end
964
965	#include "1DUTILS.h"
966	#include "1Dconv.h"
967
968

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

Download in other formats:

Original Format