Context Navigation

lmdz1d.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 34.7 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 = 59 ==> forcing_sandu = .true.
269	! initial profiles from sanduref file: see prof.inp.001
270	! SST varying with time and divergence constante: see ifa_sanduref.txt file
271	! Radiation has to be computed interactively
272	!forcing_type = 60 ==> forcing_astex = .true.
273	! initial profiles from file: see prof.inp.001
274	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
275	! Radiation has to be computed interactively
276	!forcing_type = 61 ==> forcing_armcu = .true.
277	! initial profiles from file: see prof.inp.001
278	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
279	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
280	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
281	! Radiation to be switched off
282	!
283	if (forcing_type .eq.0) THEN
284	forcing_les = .true.
285	elseif (forcing_type .eq.1) THEN
286	forcing_radconv = .true.
287	elseif (forcing_type .eq.2) THEN
288	forcing_toga = .true.
289	elseif (forcing_type .eq.3) THEN
290	forcing_GCM2SCM = .true.
291	elseif (forcing_type .eq.4) THEN
292	forcing_twpice = .true.
293	elseif (forcing_type .eq.5) THEN
294	forcing_rico = .true.
295	elseif (forcing_type .eq.6) THEN
296	forcing_amma = .true.
297	elseif (forcing_type .eq.40) THEN
298	forcing_GCSSold = .true.
299	elseif (forcing_type .eq.59) THEN
300	forcing_sandu = .true.
301	elseif (forcing_type .eq.60) THEN
302	forcing_astex = .true.
303	elseif (forcing_type .eq.61) THEN
304	forcing_armcu = .true.
305	IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
306	else
307	write (,) 'ERROR : unknown forcing_type ', forcing_type
308	stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
309	ENDIF
310	print*,"forcing type=",forcing_type
311
312	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
313	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
314	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
315	! through the common sst_forcing.
316
317	type_ts_forcing = 0
318	if (forcing_toga .or. forcing_sandu .or. forcing_astex)
319	: type_ts_forcing = 1
320
321	!---------------------------------------------------------------------
322	! Definition of the run
323	!---------------------------------------------------------------------
324
325	call conf_gcm( 99, .TRUE. , clesphy0 )
326	c-----------------------------------------------------------------------
327	c Choix du calendrier
328	c -------------------
329
330	c calend = 'earth_365d'
331	if (calend == 'earth_360d') then
332	call ioconf_calendar('360d')
333	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
334	else if (calend == 'earth_365d') then
335	call ioconf_calendar('noleap')
336	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
337	else if (calend == 'earth_366d') then
338	call ioconf_calendar('all_leap')
339	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
340	else if (calend == 'gregorian') then
341	call ioconf_calendar('gregorian') ! not to be used by normal users
342	write(,)'CALENDRIER CHOISI: Gregorien'
343	else
344	write (,) 'ERROR : unknown calendar ', calend
345	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
346	endif
347	c-----------------------------------------------------------------------
348	c
349	cc Date :
350	c La date est supposee donnee sous la forme [annee, numero du jour dans
351	c l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
352	c On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
353	c Le numero du jour est dans "day". L heure est traitee separement.
354	c La date complete est dans "daytime" (l'unite est le jour).
355	fnday=nday
356	c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
357	IF(forcing_type .EQ. 61) fnday=53100./86400.
358	c Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
359	IF(forcing_type .EQ. 6) fnday=64800./86400.
360	annee_ref = anneeref
361	mois = 1
362	day_ref = dayref
363	heure = 0.
364	itau_dyn = 0
365	itau_phy = 0
366	call ymds2ju(annee_ref,mois,day_ref,heure,day)
367	day_ini = day
368	day_end = day_ini + nday
369
370	IF (forcing_type .eq.2) THEN
371	! Convert the initial date of Toga-Coare to Julian day
372	call ymds2ju
373	$ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
374
375	ELSEIF (forcing_type .eq.4) THEN
376	! Convert the initial date of TWPICE to Julian day
377	call ymds2ju
378	$ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
379	$ ,day_ju_ini_twpi)
380	ELSEIF (forcing_type .eq.6) THEN
381	! Convert the initial date of AMMA to Julian day
382	call ymds2ju
383	$ (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma
384	$ ,day_ju_ini_amma)
385
386	ELSEIF (forcing_type .eq.59) THEN
387	! Convert the initial date of Sandu case to Julian day
388	call ymds2ju
389	$ (year_ini_sandu,mth_ini_sandu,day_ini_sandu,
390	$ time_ini*3600.,day_ju_ini_sandu)
391
392	ELSEIF (forcing_type .eq.60) THEN
393	! Convert the initial date of Astex case to Julian day
394	call ymds2ju
395	$ (year_ini_astex,mth_ini_astex,day_ini_astex,
396	$ time_ini*3600.,day_ju_ini_astex)
397
398	ELSEIF (forcing_type .eq.61) THEN
399
400	! Convert the initial date of Arm_cu case to Julian day
401	call ymds2ju
402	$ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
403	$ ,day_ju_ini_armcu)
404	ENDIF
405
406	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
407	! Print out the actual date of the beginning of the simulation :
408	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
409	print *,' Time of beginning : ',
410	$ year_print, month_print, day_print, sec_print
411
412	!---------------------------------------------------------------------
413	! Initialization of dimensions, geometry and initial state
414	!---------------------------------------------------------------------
415	call init_phys_lmdz(1,1,llm,1,(/1/))
416	call suphel
417	call initcomgeomphy
418	call infotrac_init
419
420	allocate(q(llm,nqtot)) ; q(:,:)=0.
421	allocate(dq(llm,nqtot))
422	allocate(dq_dyn(llm,nqtot))
423	allocate(d_q_adv(llm,nqtot))
424
425	c
426	c No ozone climatology need be read in this pre-initialization
427	c (phys_state_var_init is called again in physiq)
428	read_climoz = 0
429	c
430	call phys_state_var_init(read_climoz)
431
432	if (ngrid.ne.klon) then
433	print*,'stop in inifis'
434	print*,'Probleme de dimensions :'
435	print*,'ngrid = ',ngrid
436	print*,'klon = ',klon
437	stop
438	endif
439	!!!=====================================================================
440	!!! Feedback forcing values for Gateaux differentiation (al1)
441	!!!=====================================================================
442	!!! Surface Planck forcing bracketing call radiation
443	!! surf_Planck = 0.
444	!! surf_Conv = 0.
445	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
446	!!! a mettre dans le lmdz1d.def ou autre
447	!!
448	!!
449	qsol = qsolinp
450	qsurf = fq_sat(tsurf,psurf/100.)
451	rlat=xlat
452	rlon=xlon
453	day1= day_ini
454	time=daytime-day
455	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
456	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
457
458	!
459	!! mpl et jyg le 22/08/2012 :
460	!! pour que les cas a flux de surface imposes marchent
461	IF(.NOT.ok_flux_surf) THEN
462	fsens=-wtsurfrcpdrho(1)
463	flat=-wqsurfrlvttrho(1)
464	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
465	ENDIF
466	!! ok_flux_surf=.false.
467	!! fsens=-wtsurfrcpdrho(1)
468	!! flat=-wqsurfrlvttrho(1)
469	!!!!
470
471	! Vertical discretization and pressure levels at half and mid levels:
472
473	pa = 5e4
474	!! preff= 1.01325e5
475	preff = psurf
476	IF (ok_old_disvert) THEN
477	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
478	print *,'On utilise disvert0'
479	ELSE
480	call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
481	: scaleheight)
482	print *,'On utilise disvert'
483	c Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
484	c Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
485	ENDIF
486	sig_s=presnivs/preff
487	plev =ap+bp*psurf
488	play = 0.5*(plev(1:llm)+plev(2:llm+1))
489	ccc zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles
490
491	IF (forcing_type .eq. 59) THEN
492	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
493	write(,) '***********************'
494	do l = 1, llm
495	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
496	if (trouve_700 .and. play(l).le.70000) then
497	llm700=l
498	print *,'llm700,play=',llm700,play(l)/100.
499	trouve_700= .false.
500	endif
501	enddo
502	write(,) '***********************'
503	ENDIF
504
505	c
506	!=====================================================================
507	! EVENTUALLY, READ FORCING DATA :
508	!=====================================================================
509
510	#include "1D_read_forc_cases.h"
511
512	if (forcing_GCM2SCM) then
513	write (,) 'forcing_GCM2SCM not yet implemented'
514	stop 'in initialization'
515	endif ! forcing_GCM2SCM
516
517	print*,'mxcalc=',mxcalc
518	print*,'zlay=',zlay(mxcalc)
519	print*,'play=',play(mxcalc)
520
521	cAl1 pour SST forced, appellé depuis ocean_forced_noice
522	ts_cur = tsurf ! SST used in read_tsurf1d
523	!=====================================================================
524	! Initialisation de la physique :
525	!=====================================================================
526
527	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
528	!
529	! day_step, iphysiq lus dans gcm.def ci-dessus
530	! timestep: calcule ci-dessous from rday et day_step
531	! ngrid=1
532	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
533	! rday: defini dans suphel.F (86400.)
534	! day_ini: lu dans run.def (dayref)
535	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
536	! airefi,zcufi,zcvfi initialises au debut de ce programme
537	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
538	day_step = float(nsplit_phys)*day_step/float(iphysiq)
539	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
540	timestep =rday/day_step
541	dtime_frcg = timestep
542	!
543	zcufi=airefi
544	zcvfi=airefi
545	!
546	rlat_rad(:)=rlat(:)*rpi/180.
547	rlon_rad(:)=rlon(:)*rpi/180.
548
549	call iniphysiq(ngrid,llm,rday,day_ini,timestep,
550	. rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
551	print*,'apres iniphysiq'
552
553	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
554	co2_ppm= 330.0
555	solaire=1370.0
556
557	! Ecriture du startphy avant le premier appel a la physique.
558	! On le met juste avant pour avoir acces a tous les champs
559	! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
560
561	if (ok_writedem) then
562
563	!--------------------------------------------------------------------------
564	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
565	! need : qsol fder snow qsurf evap rugos agesno ftsoil
566	!--------------------------------------------------------------------------
567
568	type_ocean = "force"
569	run_off_lic_0(1) = restart_runoff
570	call fonte_neige_init(run_off_lic_0)
571
572	fder=0.
573	snsrf(1,:)=0. ! couverture de neige des sous surface
574	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
575	evap=0.
576	frugs(1,:)=rugos ! couverture de neige des sous surface
577	agesno = xagesno
578	tsoil(:,:,:)=tsurf
579	!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
580	! tsoil(1,1,1)=299.18
581	! tsoil(1,2,1)=300.08
582	! tsoil(1,3,1)=301.88
583	! tsoil(1,4,1)=305.48
584	! tsoil(1,5,1)=308.00
585	! tsoil(1,6,1)=308.00
586	! tsoil(1,7,1)=308.00
587	! tsoil(1,8,1)=308.00
588	! tsoil(1,9,1)=308.00
589	! tsoil(1,10,1)=308.00
590	! tsoil(1,11,1)=308.00
591	!-----------------------------------------------------------------------
592	call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,
593	& evap, frugs, agesno, tsoil)
594
595	!------------------ prepare limit conditions for limit.nc -----------------
596	!-- Ocean force
597
598	print*,'avant phyredem'
599	pctsrf(1,:)=0.
600	if (nat_surf.eq.0.) then
601	pctsrf(1,is_oce)=1.
602	pctsrf(1,is_ter)=0.
603	else
604	pctsrf(1,is_oce)=0.
605	pctsrf(1,is_ter)=1.
606	end if
607
608	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf
609	$ ,pctsrf(1,is_oce),pctsrf(1,is_ter)
610
611	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
612	zpic = zpicinp
613	ftsol=tsurf
614	falb1 = albedo
615	falb2 = albedo
616	rugoro=rugos
617	t_ancien(1,:)=temp(:)
618	q_ancien(1,:)=q(:,1)
619	pbl_tke=1.e-8
620
621	rain_fall=0.
622	snow_fall=0.
623	solsw=0.
624	sollw=0.
625	radsol=0.
626	rnebcon=0.
627	ratqs=0.
628	clwcon=0.
629	zmea=0.
630	zstd=0.
631	zsig=0.
632	zgam=0.
633	zval=0.
634	zthe=0.
635	sig1=0.
636	w01=0.
637	u_ancien(1,:)=u(:)
638	v_ancien(1,:)=v(:)
639
640	!------------------------------------------------------------------------
641	! Make file containing restart for the physics (startphy.nc)
642	!
643	! NB: List of the variables to be written by phyredem (via put_field):
644	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
645	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
646	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
647	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
648	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
649	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
650	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,sig1,w01
651	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
652	!------------------------------------------------------------------------
653	CAl1 =============== restart option ==========================
654	if (.not.restart) then
655	call phyredem ("startphy.nc")
656	else
657	c (desallocations)
658	print*,'callin surf final'
659	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,
660	& evap, frugs, agesno, tsoil)
661	print*,'after surf final'
662	CALL fonte_neige_final(run_off_lic_0)
663	endif
664
665	ok_writedem=.false.
666	print*,'apres phyredem'
667
668	endif ! ok_writedem
669
670	!------------------------------------------------------------------------
671	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
672	! --------------------------------------------------
673	! NB: List of the variables to be written in limit.nc
674	! (by writelim.F, subroutine of 1DUTILS.h):
675	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
676	! phy_fter,phy_foce,phy_flic,phy_fsic)
677	!------------------------------------------------------------------------
678	do i=1,yd
679	phy_nat(i) = nat_surf
680	phy_alb(i) = albedo
681	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
682	phy_rug(i) = rugos
683	phy_fter(i) = pctsrf(1,is_ter)
684	phy_foce(i) = pctsrf(1,is_oce)
685	phy_fsic(i) = pctsrf(1,is_sic)
686	phy_flic(i) = pctsrf(1,is_lic)
687	enddo
688
689	C fabrication de limit.nc
690	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,
691	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
692
693
694	call phys_state_var_end
695	cAl1
696	if (restart) then
697	print*,'call to restart dyn 1d'
698	Call dyn1deta0("start1dyn.nc",
699	& plev,play,phi,phis,presnivs,
700	& u,v,temp,q,omega2)
701
702	print*,'fnday,annee_ref,day_ref,day_ini',
703	& fnday,annee_ref,day_ref,day_ini
704	c** call ymds2ju(annee_ref,mois,day_ini,heure,day)
705	day = day_ini
706	day_end = day_ini + nday
707	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
708
709	! Print out the actual date of the beginning of the simulation :
710	call ju2ymds(daytime, an, mois, jour, heure)
711	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
712
713	day = int(daytime)
714	time=daytime-day
715
716	print,'*** intialised fields from restart1dyn *****'
717	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
718	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
719	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
720	c raz for safety
721	do l=1,llm
722	dq_dyn(l,1) = 0.
723	enddo
724	endif
725	!Al1 ================ end restart =================================
726	IF (ecrit_slab_oc.eq.1) then
727	open(97,file='div_slab.dat',STATUS='UNKNOWN')
728	elseif (ecrit_slab_oc.eq.0) then
729	open(97,file='div_slab.dat',STATUS='OLD')
730	endif
731	!=====================================================================
732	! START OF THE TEMPORAL LOOP :
733	!=====================================================================
734
735	do while(it.le.nint(fnday*day_step))
736
737	if (prt_level.ge.1) then
738	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',
739	. it,day,time,nint(fnday*day_step)
740	print*,'PAS DE TEMPS ',timestep
741	endif
742	!Al1 demande de restartphy.nc
743	if (it.eq.nint(fnday*day_step)) lastcall=.True.
744
745	!---------------------------------------------------------------------
746	! Interpolation of forcings in time and onto model levels
747	!---------------------------------------------------------------------
748
749	#include "1D_interp_cases.h"
750
751	if (forcing_GCM2SCM) then
752	write (,) 'forcing_GCM2SCM not yet implemented'
753	stop 'in time loop'
754	endif ! forcing_GCM2SCM
755
756	!---------------------------------------------------------------------
757	! Geopotential :
758	!---------------------------------------------------------------------
759
760	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
761	do l = 1, llm-1
762	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1))
763	. (play(l)-play(l+1))/(play(l)+play(l+1))
764	enddo
765
766	!---------------------------------------------------------------------
767	! Listing output for debug prt_level>=1
768	!---------------------------------------------------------------------
769	if (prt_level>=1) then
770	print *,' avant physiq : -------- day time ',day,time
771	write(,) 'firstcall,lastcall,phis',
772	: firstcall,lastcall,phis
773	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',
774	: 'presniv','plev','play','phi'
775	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,
776	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
777	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',
778	: 'presniv','u','v','temp','q1','q2','omega2'
779	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,
780	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
781	endif
782
783	!---------------------------------------------------------------------
784	! Call physiq :
785	!---------------------------------------------------------------------
786
787	call physiq(ngrid,llm,
788	: firstcall,lastcall,
789	: day,time,timestep,
790	: plev,play,phi,phis,presnivs,clesphy0,
791	: u,v,temp,q,omega2,
792	: du_phys,dv_phys,dt_phys,dq,dpsrf,
793	: dudyn,PVteta)
794	firstcall=.false.
795
796	!---------------------------------------------------------------------
797	! Listing output for debug prt_level>=1
798	!---------------------------------------------------------------------
799	if (prt_level>=1) then
800	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',
801	: 'presniv','plev','play','phi'
802	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,
803	: presnivs(l),plev(l),play(l),phi(l),l=1,llm)
804	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',
805	: 'presniv','u','v','temp','q1','q2','omega2'
806	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,
807	: presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
808	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',
809	: 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
810	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,
811	: presnivs(l),86400du_phys(l),86400dv_phys(l),
812	: 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
813	write(,) 'dpsrf',dpsrf
814	endif
815	!---------------------------------------------------------------------
816	! Add physical tendencies :
817	!---------------------------------------------------------------------
818
819	fcoriolis=2.sin(rpixlat/180.)*romega
820	if (forcing_radconv) then
821	fcoriolis=0.0
822	dt_cooling=0.0
823	d_th_adv=0.0
824	d_q_adv=0.0
825	endif
826
827	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice
828	: .or.forcing_amma) then
829	fcoriolis=0.0 ; ug=0. ; vg=0.
830	endif
831	if(forcing_rico) then
832	dt_cooling=0.
833	endif
834
835	print*, 'fcoriolis ', fcoriolis, xlat
836
837	du_age(1:mxcalc)=
838	: fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
839	dv_age(1:mxcalc)=
840	: -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
841
842	!!!!!!!!!!!!!!!!!!!!!!!!
843	! Geostrophic wind
844	!!!!!!!!!!!!!!!!!!!!!!!!
845	sfdt = sin(0.5fcoriolistimestep)
846	cfdt = cos(0.5fcoriolistimestep)
847	!
848	du_age(1:mxcalc)= -2.sfdt/timestep
849	: (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -
850	: cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
851	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
852	!
853	dv_age(1:mxcalc)= -2.sfdt/timestep
854	: (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +
855	: sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
856	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
857	!
858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859	! call writefield_phy('dv_age' ,dv_age,llm)
860	! call writefield_phy('du_age' ,du_age,llm)
861	! call writefield_phy('du_phys' ,du_phys,llm)
862	! call writefield_phy('u_tend' ,u,llm)
863	! call writefield_phy('u_g' ,ug,llm)
864	!
865	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
866	!! Increment state variables
867	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
868	u(1:mxcalc)=u(1:mxcalc) + timestep*(
869	: du_phys(1:mxcalc)
870	: +du_age(1:mxcalc) )
871	v(1:mxcalc)=v(1:mxcalc) + timestep*(
872	: dv_phys(1:mxcalc)
873	: +dv_age(1:mxcalc) )
874	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
875	: dq(1:mxcalc,:)
876	: +d_q_adv(1:mxcalc,:) )
877
878	if (prt_level.ge.1) then
879	print *,
880	: 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',
881	: temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
882	print*,du_phys
883	print*, v
884	print*, vg
885	endif
886
887	temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
888	. dt_phys(1:mxcalc)
889	. +d_th_adv(1:mxcalc)
890	. +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
891
892	teta=temp(pzero/play)*rkappa
893	!
894	!---------------------------------------------------------------------
895	! Nudge soil temperature if requested
896	!---------------------------------------------------------------------
897
898	IF (nudge_tsoil) THEN
899	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)
900	. -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
901	ENDIF
902
903	!---------------------------------------------------------------------
904	! Add large-scale tendencies (advection, etc) :
905	!---------------------------------------------------------------------
906
907	ccc nrlmd
908	ccc tmpvar=teta
909	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
910	ccc
911	ccc teta(1:mxcalc)=tmpvar(1:mxcalc)
912	ccc tmpvar(:)=q(:,1)
913	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
914	ccc q(1:mxcalc,1)=tmpvar(1:mxcalc)
915	ccc tmpvar(:)=q(:,2)
916	ccc call advect_vert(llm,omega,timestep,tmpvar,plev)
917	ccc q(1:mxcalc,2)=tmpvar(1:mxcalc)
918
919	!---------------------------------------------------------------------
920	! Air temperature :
921	!---------------------------------------------------------------------
922	if (lastcall) then
923	print*,'Pas de temps final ',it
924	call ju2ymds(daytime, an, mois, jour, heure)
925	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
926	endif
927
928	c incremente day time
929	c print*,'daytime bef',daytime,1./day_step
930	daytime = daytime+1./day_step
931	!Al1dbg
932	day = int(daytime+0.1/day_step)
933	c time = max(daytime-day,0.0)
934	cAl1&jyg: correction de bug
935	ccc time = real(mod(it,day_step))/day_step
936	time = time_ini/24.+real(mod(it,day_step))/day_step
937	c print*,'daytime nxt time',daytime,time
938	it=it+1
939
940	enddo
941
942	!Al1
943	if (ecrit_slab_oc.ne.-1) close(97)
944
945	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
946	! -------------------------------------
947	call dyn1dredem("restart1dyn.nc",
948	: plev,play,phi,phis,presnivs,
949	: u,v,temp,q,omega2)
950
951	CALL abort_gcm ('lmdz1d ','The End ',0)
952
953	end
954
955	#include "1DUTILS.h"
956	#include "1Dconv.h"
957
958

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

Download in other formats:

Original Format