Context Navigation

lmdz1d.F90 @ 2119

Last change on this file since 2119 was 2117, checked in by fhourdin, 10 years ago

Modification de la lecture du cas 1D AMMA en vue d'une generatlisation
de la specification des forcages decidee dans le cadre de DEPHY.
Passage a une allocation dynamique des fichiers pour permettre
de lire la dimension des variables dans le fichier de forcage.
Marche avec AMMA et le noveau cas IHOP.

Change of 1D cases forcing. Introduction of allocated variable.
For AMMA and IHOP.

File size: 36.9 KB

Line
1	#ifdef CPP_1D
2	#include "../dyn3d/mod_const_mpi.F90"
3	#include "../dyn3d_common/control_mod.F90"
4	#include "../dyn3d_common/infotrac.F90"
5	#include "../dyn3d_common/disvert.F90"
6
7
8	PROGRAM lmdz1d
9
10	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
11	use phys_state_var_mod
12	use comgeomphy
13	use dimphy
14	use surface_data, only : type_ocean,ok_veget
15	use pbl_surface_mod, only : ftsoil, pbl_surface_init, &
16	& pbl_surface_final
17	use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
18
19	use infotrac ! new
20	use control_mod
21	USE indice_sol_mod
22	USE phyaqua_mod
23	USE mod_1D_cases_read
24
25	implicit none
26	#include "dimensions.h"
27	#include "YOMCST.h"
28	#include "temps.h"
29	!!#include "control.h"
30	#include "iniprint.h"
31	#include "clesphys.h"
32	#include "dimsoil.h"
33	!#include "indicesol.h"
34
35	#include "comvert.h"
36	#include "compar1d.h"
37	#include "flux_arp.h"
38	#include "tsoilnudge.h"
39	#include "fcg_gcssold.h"
40	!!!#include "fbforcing.h"
41
42	!=====================================================================
43	! DECLARATIONS
44	!=====================================================================
45
46	!---------------------------------------------------------------------
47	! Externals
48	!---------------------------------------------------------------------
49	external fq_sat
50	real fq_sat
51
52	!---------------------------------------------------------------------
53	! Arguments d' initialisations de la physique (USER DEFINE)
54	!---------------------------------------------------------------------
55
56	integer, parameter :: ngrid=1
57	real :: zcufi = 1.
58	real :: zcvfi = 1.
59
60	!- real :: nat_surf
61	!- logical :: ok_flux_surf
62	!- real :: fsens
63	!- real :: flat
64	!- real :: tsurf
65	!- real :: rugos
66	!- real :: qsol(1:2)
67	!- real :: qsurf
68	!- real :: psurf
69	!- real :: zsurf
70	!- real :: albedo
71	!-
72	!- real :: time = 0.
73	!- real :: time_ini
74	!- real :: xlat
75	!- real :: xlon
76	!- real :: wtsurf
77	!- real :: wqsurf
78	!- real :: restart_runoff
79	!- real :: xagesno
80	!- real :: qsolinp
81	!- real :: zpicinp
82	!-
83	real :: fnday
84	real :: day, daytime
85	real :: day1
86	real :: heure
87	integer :: jour
88	integer :: mois
89	integer :: an
90
91	!---------------------------------------------------------------------
92	! Declarations related to forcing and initial profiles
93	!---------------------------------------------------------------------
94
95	integer :: kmax = llm
96	integer llm700,nq1,nq2
97	INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
98	real timestep, frac
99	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
100	real uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
101	real ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
102	real dqtdxls(nlev_max),dqtdyls(nlev_max)
103	real dqtdtls(nlev_max),thlpcar(nlev_max)
104	real qprof(nlev_max,nqmx)
105
106	! integer :: forcing_type
107	logical :: forcing_les = .false.
108	logical :: forcing_armcu = .false.
109	logical :: forcing_rico = .false.
110	logical :: forcing_radconv = .false.
111	logical :: forcing_toga = .false.
112	logical :: forcing_twpice = .false.
113	logical :: forcing_amma = .false.
114	logical :: forcing_GCM2SCM = .false.
115	logical :: forcing_GCSSold = .false.
116	logical :: forcing_sandu = .false.
117	logical :: forcing_astex = .false.
118	logical :: forcing_fire = .false.
119	integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
120	! (cf read_tsurf1d.F)
121
122	!vertical advection computation
123	! real d_t_z(llm), d_q_z(llm)
124	! real d_t_dyn_z(llm), d_q_dyn_z(llm)
125	! real zz(llm)
126	! real zfact
127
128	!flag forcings
129	logical :: nudge_wind=.true.
130	logical :: nudge_thermo=.false.
131	logical :: cptadvw=.true.
132	!=====================================================================
133	! DECLARATIONS FOR EACH CASE
134	!=====================================================================
135	!
136	#include "1D_decl_cases.h"
137	!
138	!---------------------------------------------------------------------
139	! Declarations related to vertical discretization:
140	!---------------------------------------------------------------------
141	real :: pzero=1.e5
142	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
143	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
144
145	!---------------------------------------------------------------------
146	! Declarations related to variables
147	!---------------------------------------------------------------------
148
149	real :: phi(llm)
150	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
151	real :: rlat_rad(1),rlon_rad(1)
152	real :: omega(llm+1),omega2(llm),rho(llm+1)
153	real :: ug(llm),vg(llm),fcoriolis
154	real :: sfdt, cfdt
155	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
156	real :: dt_dyn(llm)
157	real :: dt_cooling(llm),d_th_adv(llm)
158	real :: alpha
159	real :: ttt
160
161	REAL, ALLOCATABLE, DIMENSION(:,:):: q
162	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
163	REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
164	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
165	! REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
166
167	!---------------------------------------------------------------------
168	! Initialization of surface variables
169	!---------------------------------------------------------------------
170	real :: run_off_lic_0(1)
171	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
172	real :: evap(1,nbsrf),frugs(1,nbsrf)
173	real :: tsoil(1,nsoilmx,nbsrf)
174	real :: agesno(1,nbsrf)
175
176	!---------------------------------------------------------------------
177	! Call to phyredem
178	!---------------------------------------------------------------------
179	logical :: ok_writedem =.true.
180
181	!---------------------------------------------------------------------
182	! Call to physiq
183	!---------------------------------------------------------------------
184	integer, parameter :: longcles=20
185	logical :: firstcall=.true.
186	logical :: lastcall=.false.
187	real :: phis = 0.0
188	real :: clesphy0(longcles) = 0.0
189	real :: dpsrf
190
191	!---------------------------------------------------------------------
192	! Initializations of boundary conditions
193	!---------------------------------------------------------------------
194	integer, parameter :: yd = 360
195	real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
196	real :: phy_alb (yd) ! Albedo land only (old value condsurf_jyg=0.3)
197	real :: phy_sst (yd) ! SST (will not be used; cf read_tsurf1d.F)
198	real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
199	real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
200	real :: phy_ice (yd) = 0.0 ! Fraction de glace
201	real :: phy_fter(yd) = 0.0 ! Fraction de terre
202	real :: phy_foce(yd) = 0.0 ! Fraction de ocean
203	real :: phy_fsic(yd) = 0.0 ! Fraction de glace
204	real :: phy_flic(yd) = 0.0 ! Fraction de glace
205
206	!---------------------------------------------------------------------
207	! Fichiers et d'autres variables
208	!---------------------------------------------------------------------
209	integer :: k,l,i,it=1,mxcalc
210	integer jjmp1
211	parameter (jjmp1=jjm+1-1/jjm)
212	REAL dudyn(iim+1,jjmp1,llm)
213	INTEGER read_climoz
214	!Al1
215	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
216	data ecrit_slab_oc/-1/
217
218	!=====================================================================
219	! INITIALIZATIONS
220	!=====================================================================
221	! Initialization of Common turb_forcing
222	dtime_frcg = 0.
223	Turb_fcg_gcssold=.false.
224	hthturb_gcssold = 0.
225	hqturb_gcssold = 0.
226
227	!---------------------------------------------------------------------
228	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
229	!---------------------------------------------------------------------
230	!Al1
231	call conf_unicol
232	!Al1 moves this gcssold var from common fcg_gcssold to
233	Turb_fcg_gcssold = xTurb_fcg_gcssold
234	! --------------------------------------------------------------------
235	close(1)
236	!Al1
237	write(,) 'lmdz1d.def lu => unicol.def'
238
239	! forcing_type defines the way the SCM is forced:
240	!forcing_type = 0 ==> forcing_les = .true.
241	! initial profiles from file prof.inp.001
242	! no forcing by LS convergence ;
243	! surface temperature imposed ;
244	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
245	!forcing_type = 1 ==> forcing_radconv = .true.
246	! idem forcing_type = 0, but the imposed radiative cooling
247	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
248	! then there is no radiative cooling at all)
249	!forcing_type = 2 ==> forcing_toga = .true.
250	! initial profiles from TOGA-COARE IFA files
251	! LS convergence and SST imposed from TOGA-COARE IFA files
252	!forcing_type = 3 ==> forcing_GCM2SCM = .true.
253	! initial profiles from the GCM output
254	! LS convergence imposed from the GCM output
255	!forcing_type = 4 ==> forcing_twpice = .true.
256	! initial profiles from TWP-ICE cdf file
257	! LS convergence, omega and SST imposed from TWP-ICE files
258	!forcing_type = 5 ==> forcing_rico = .true.
259	! initial profiles from RICO files
260	! LS convergence imposed from RICO files
261	!forcing_type = 6 ==> forcing_amma = .true.
262	! initial profiles from AMMA nc file
263	! LS convergence, omega and surface fluxes imposed from AMMA file
264	!forcing_type = 40 ==> forcing_GCSSold = .true.
265	! initial profile from GCSS file
266	! LS convergence imposed from GCSS file
267	!forcing_type = 50 ==> forcing_fire = .true.
268	! forcing from fire.nc
269	!forcing_type = 59 ==> forcing_sandu = .true.
270	! initial profiles from sanduref file: see prof.inp.001
271	! SST varying with time and divergence constante: see ifa_sanduref.txt file
272	! Radiation has to be computed interactively
273	!forcing_type = 60 ==> forcing_astex = .true.
274	! initial profiles from file: see prof.inp.001
275	! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
276	! Radiation has to be computed interactively
277	!forcing_type = 61 ==> forcing_armcu = .true.
278	! initial profiles from file: see prof.inp.001
279	! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
280	! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
281	! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
282	! Radiation to be switched off
283	!
284	if (forcing_type .eq.0) THEN
285	forcing_les = .true.
286	elseif (forcing_type .eq.1) THEN
287	forcing_radconv = .true.
288	elseif (forcing_type .eq.2) THEN
289	forcing_toga = .true.
290	elseif (forcing_type .eq.3) THEN
291	forcing_GCM2SCM = .true.
292	elseif (forcing_type .eq.4) THEN
293	forcing_twpice = .true.
294	elseif (forcing_type .eq.5) THEN
295	forcing_rico = .true.
296	elseif (forcing_type .eq.6) THEN
297	forcing_amma = .true.
298	elseif (forcing_type .eq.40) THEN
299	forcing_GCSSold = .true.
300	elseif (forcing_type .eq.50) THEN
301	forcing_fire = .true.
302	elseif (forcing_type .eq.59) THEN
303	forcing_sandu = .true.
304	elseif (forcing_type .eq.60) THEN
305	forcing_astex = .true.
306	elseif (forcing_type .eq.61) THEN
307	forcing_armcu = .true.
308	IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
309	else
310	write (,) 'ERROR : unknown forcing_type ', forcing_type
311	stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
312	ENDIF
313	print*,"forcing type=",forcing_type
314
315	! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
316	! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
317	! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
318	! through the common sst_forcing.
319
320	type_ts_forcing = 0
321	if (forcing_toga.or.forcing_sandu.or.forcing_astex) &
322	& type_ts_forcing = 1
323
324	!---------------------------------------------------------------------
325	! Definition of the run
326	!---------------------------------------------------------------------
327
328	call conf_gcm( 99, .TRUE. , clesphy0 )
329	!-----------------------------------------------------------------------
330	! Choix du calendrier
331	! -------------------
332
333	! calend = 'earth_365d'
334	if (calend == 'earth_360d') then
335	call ioconf_calendar('360d')
336	write(,)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
337	else if (calend == 'earth_365d') then
338	call ioconf_calendar('noleap')
339	write(,)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
340	else if (calend == 'earth_366d') then
341	call ioconf_calendar('all_leap')
342	write(,)'CALENDRIER CHOISI: Terrestre bissextile'
343	else if (calend == 'gregorian') then
344	call ioconf_calendar('gregorian') ! not to be used by normal users
345	write(,)'CALENDRIER CHOISI: Gregorien'
346	else
347	write (,) 'ERROR : unknown calendar ', calend
348	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
349	endif
350	!-----------------------------------------------------------------------
351	!
352	!c Date :
353	! La date est supposee donnee sous la forme [annee, numero du jour dans
354	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
355	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
356	! Le numero du jour est dans "day". L heure est traitee separement.
357	! La date complete est dans "daytime" (l'unite est le jour).
358	if (nday>0) then
359	fnday=nday
360	else
361	fnday=-nday/float(day_step)
362	endif
363
364	! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
365	IF(forcing_type .EQ. 61) fnday=53100./86400.
366	! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
367	IF(forcing_type .EQ. 6) fnday=64800./86400.
368	annee_ref = anneeref
369	mois = 1
370	day_ref = dayref
371	heure = 0.
372	itau_dyn = 0
373	itau_phy = 0
374	call ymds2ju(annee_ref,mois,day_ref,heure,day)
375	day_ini = int(day)
376	day_end = day_ini + fnday
377
378	IF (forcing_type .eq.2) THEN
379	! Convert the initial date of Toga-Coare to Julian day
380	call ymds2ju &
381	& (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
382
383	ELSEIF (forcing_type .eq.4) THEN
384	! Convert the initial date of TWPICE to Julian day
385	call ymds2ju &
386	& (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi &
387	& ,day_ju_ini_twpi)
388	ELSEIF (forcing_type .eq.6) THEN
389	! Convert the initial date of AMMA to Julian day
390	call ymds2ju &
391	& (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma &
392	& ,day_ju_ini_amma)
393
394	ELSEIF (forcing_type .eq.59) THEN
395	! Convert the initial date of Sandu case to Julian day
396	call ymds2ju &
397	& (year_ini_sandu,mth_ini_sandu,day_ini_sandu, &
398	& time_ini*3600.,day_ju_ini_sandu)
399
400	ELSEIF (forcing_type .eq.60) THEN
401	! Convert the initial date of Astex case to Julian day
402	call ymds2ju &
403	& (year_ini_astex,mth_ini_astex,day_ini_astex, &
404	& time_ini*3600.,day_ju_ini_astex)
405
406	ELSEIF (forcing_type .eq.61) THEN
407
408	! Convert the initial date of Arm_cu case to Julian day
409	call ymds2ju &
410	& (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu &
411	& ,day_ju_ini_armcu)
412	ENDIF
413
414	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
415	! Print out the actual date of the beginning of the simulation :
416	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
417	print *,' Time of beginning : ', &
418	& year_print, month_print, day_print, sec_print
419
420	!---------------------------------------------------------------------
421	! Initialization of dimensions, geometry and initial state
422	!---------------------------------------------------------------------
423	call init_phys_lmdz(1,1,llm,1,(/1/))
424	call suphel
425	call initcomgeomphy
426	call infotrac_init
427
428	if (nqtot>nqmx) STOP'Augmenter nqmx dans lmdz1d.F'
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	! allocate(d_th_adv(llm))
434
435	!
436	! No ozone climatology need be read in this pre-initialization
437	! (phys_state_var_init is called again in physiq)
438	read_climoz = 0
439	!
440	call phys_state_var_init(read_climoz)
441
442	if (ngrid.ne.klon) then
443	print*,'stop in inifis'
444	print*,'Probleme de dimensions :'
445	print*,'ngrid = ',ngrid
446	print*,'klon = ',klon
447	stop
448	endif
449	!!!=====================================================================
450	!!! Feedback forcing values for Gateaux differentiation (al1)
451	!!!=====================================================================
452	!!! Surface Planck forcing bracketing call radiation
453	!! surf_Planck = 0.
454	!! surf_Conv = 0.
455	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
456	!!! a mettre dans le lmdz1d.def ou autre
457	!!
458	!!
459	qsol = qsolinp
460	qsurf = fq_sat(tsurf,psurf/100.)
461	rlat=xlat
462	rlon=xlon
463	day1= day_ini
464	time=daytime-day
465	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
466	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
467
468	!
469	!! mpl et jyg le 22/08/2012 :
470	!! pour que les cas a flux de surface imposes marchent
471	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
472	fsens=-wtsurfrcpdrho(1)
473	flat=-wqsurfrlvttrho(1)
474	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
475	ENDIF
476	print*,'Flux sol ',fsens,flat
477	!! ok_flux_surf=.false.
478	!! fsens=-wtsurfrcpdrho(1)
479	!! flat=-wqsurfrlvttrho(1)
480	!!!!
481
482	! Vertical discretization and pressure levels at half and mid levels:
483
484	pa = 5e4
485	!! preff= 1.01325e5
486	preff = psurf
487	IF (ok_old_disvert) THEN
488	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
489	print *,'On utilise disvert0'
490	ELSE
491	call disvert()
492	print *,'On utilise disvert'
493	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
494	! 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	!cc 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	!
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	!Al1 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	!Al1 =============== restart option ==========================
664	if (.not.restart) then
665	call phyredem ("startphy.nc")
666	else
667	! (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	! 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	!Al1
706	if (restart) then
707	print*,'call to restart dyn 1d'
708	Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs, &
709	& u,v,temp,q,omega2)
710
711	print*,'fnday,annee_ref,day_ref,day_ini', &
712	& fnday,annee_ref,day_ref,day_ini
713	!** call ymds2ju(annee_ref,mois,day_ini,heure,day)
714	day = day_ini
715	day_end = day_ini + nday
716	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
717
718	! Print out the actual date of the beginning of the simulation :
719	call ju2ymds(daytime, an, mois, jour, heure)
720	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
721
722	day = int(daytime)
723	time=daytime-day
724
725	print,'*** intialised fields from restart1dyn *****'
726	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
727	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
728	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
729	! raz for safety
730	do l=1,llm
731	dq_dyn(l,1) = 0.
732	enddo
733	endif
734	!Al1 ================ end restart =================================
735	IF (ecrit_slab_oc.eq.1) then
736	open(97,file='div_slab.dat',STATUS='UNKNOWN')
737	elseif (ecrit_slab_oc.eq.0) then
738	open(97,file='div_slab.dat',STATUS='OLD')
739	endif
740	!=====================================================================
741	! START OF THE TEMPORAL LOOP :
742	!=====================================================================
743
744	do while(it.le.nint(fnday*day_step))
745
746	if (prt_level.ge.1) then
747	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
748	& it,day,time,nint(fnday*day_step)
749	print*,'PAS DE TEMPS ',timestep
750	endif
751	!Al1 demande de restartphy.nc
752	if (it.eq.nint(fnday*day_step)) lastcall=.True.
753
754	!---------------------------------------------------------------------
755	! Interpolation of forcings in time and onto model levels
756	!---------------------------------------------------------------------
757
758	#include "1D_interp_cases.h"
759
760	if (forcing_GCM2SCM) then
761	write (,) 'forcing_GCM2SCM not yet implemented'
762	stop 'in time loop'
763	endif ! forcing_GCM2SCM
764
765	!---------------------------------------------------------------------
766	! Geopotential :
767	!---------------------------------------------------------------------
768
769	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
770	do l = 1, llm-1
771	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
772	& (play(l)-play(l+1))/(play(l)+play(l+1))
773	enddo
774
775	!---------------------------------------------------------------------
776	! Listing output for debug prt_level>=1
777	!---------------------------------------------------------------------
778	if (prt_level>=1) then
779	print *,' avant physiq : -------- day time ',day,time
780	write(,) 'firstcall,lastcall,phis', &
781	& firstcall,lastcall,phis
782	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l', &
783	& 'presniv','plev','play','phi'
784	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l, &
785	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
786	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l', &
787	& 'presniv','u','v','temp','q1','q2','omega2'
788	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l, &
789	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
790	endif
791
792	!---------------------------------------------------------------------
793	! Call physiq :
794	!---------------------------------------------------------------------
795
796	call physiq(ngrid,llm, &
797	& firstcall,lastcall, &
798	& day,time,timestep, &
799	& plev,play,phi,phis,presnivs,clesphy0, &
800	& u,v,temp,q,omega2, &
801	& du_phys,dv_phys,dt_phys,dq,dpsrf, &
802	& dudyn)
803	firstcall=.false.
804
805	!---------------------------------------------------------------------
806	! Listing output for debug prt_level>=1
807	!---------------------------------------------------------------------
808	if (prt_level>=1) then
809	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l', &
810	& 'presniv','plev','play','phi'
811	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l, &
812	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
813	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l', &
814	& 'presniv','u','v','temp','q1','q2','omega2'
815	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l, &
816	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
817	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l', &
818	& 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
819	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l, &
820	& presnivs(l),86400du_phys(l),86400dv_phys(l), &
821	& 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
822	write(,) 'dpsrf',dpsrf
823	endif
824	!---------------------------------------------------------------------
825	! Add physical tendencies :
826	!---------------------------------------------------------------------
827
828	fcoriolis=2.sin(rpixlat/180.)*romega
829	if (forcing_radconv .or. forcing_fire) then
830	fcoriolis=0.0
831	dt_cooling=0.0
832	d_th_adv=0.0
833	d_q_adv=0.0
834	endif
835
836	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice &
837	& .or.forcing_amma) then
838	fcoriolis=0.0 ; ug=0. ; vg=0.
839	endif
840	if(forcing_rico) then
841	dt_cooling=0.
842	endif
843
844	print*, 'fcoriolis ', fcoriolis, xlat
845
846	du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
847	dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
848
849	!!!!!!!!!!!!!!!!!!!!!!!!
850	! Geostrophic wind
851	!!!!!!!!!!!!!!!!!!!!!!!!
852	sfdt = sin(0.5fcoriolistimestep)
853	cfdt = cos(0.5fcoriolistimestep)
854	!
855	du_age(1:mxcalc)= -2.sfdt/timestep &
856	& (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - &
857	& cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
858	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
859	!
860	dv_age(1:mxcalc)= -2.sfdt/timestep &
861	& (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + &
862	& sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
863	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
864	!
865	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
866	! call writefield_phy('dv_age' ,dv_age,llm)
867	! call writefield_phy('du_age' ,du_age,llm)
868	! call writefield_phy('du_phys' ,du_phys,llm)
869	! call writefield_phy('u_tend' ,u,llm)
870	! call writefield_phy('u_g' ,ug,llm)
871	!
872	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
873	!! Increment state variables
874	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
875	! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
876	! au dessus de 700hpa, on relaxe vers les profils initiaux
877	if (forcing_sandu .OR. forcing_astex) then
878	#include "1D_nudge_sandu_astex.h"
879	else
880	u(1:mxcalc)=u(1:mxcalc) + timestep*( &
881	& du_phys(1:mxcalc) &
882	& +du_age(1:mxcalc) )
883	v(1:mxcalc)=v(1:mxcalc) + timestep*( &
884	& dv_phys(1:mxcalc) &
885	& +dv_age(1:mxcalc) )
886	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( &
887	& dq(1:mxcalc,:) &
888	& +d_q_adv(1:mxcalc,:) )
889
890	if (prt_level.ge.1) then
891	print *, &
892	& 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ', &
893	& temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
894	print*,du_phys
895	print*, v
896	print*, vg
897	endif
898
899	temp(1:mxcalc)=temp(1:mxcalc)+timestep*( &
900	& dt_phys(1:mxcalc) &
901	& +d_th_adv(1:mxcalc) &
902	& +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
903
904	endif ! forcing_sandu or forcing_astex
905
906	teta=temp(pzero/play)*rkappa
907	!
908	!---------------------------------------------------------------------
909	! Nudge soil temperature if requested
910	!---------------------------------------------------------------------
911
912	IF (nudge_tsoil) THEN
913	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) &
914	& -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
915	ENDIF
916
917	!---------------------------------------------------------------------
918	! Add large-scale tendencies (advection, etc) :
919	!---------------------------------------------------------------------
920
921	!cc nrlmd
922	!cc tmpvar=teta
923	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
924	!cc
925	!cc teta(1:mxcalc)=tmpvar(1:mxcalc)
926	!cc tmpvar(:)=q(:,1)
927	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
928	!cc q(1:mxcalc,1)=tmpvar(1:mxcalc)
929	!cc tmpvar(:)=q(:,2)
930	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
931	!cc q(1:mxcalc,2)=tmpvar(1:mxcalc)
932
933	!---------------------------------------------------------------------
934	! Air temperature :
935	!---------------------------------------------------------------------
936	if (lastcall) then
937	print*,'Pas de temps final ',it
938	call ju2ymds(daytime, an, mois, jour, heure)
939	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
940	endif
941
942	! incremente day time
943	! print*,'daytime bef',daytime,1./day_step
944	daytime = daytime+1./day_step
945	!Al1dbg
946	day = int(daytime+0.1/day_step)
947	! time = max(daytime-day,0.0)
948	!Al1&jyg: correction de bug
949	!cc time = real(mod(it,day_step))/day_step
950	time = time_ini/24.+real(mod(it,day_step))/day_step
951	! print*,'daytime nxt time',daytime,time
952	it=it+1
953
954	enddo
955
956	!Al1
957	if (ecrit_slab_oc.ne.-1) close(97)
958
959	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
960	! -------------------------------------
961	call dyn1dredem("restart1dyn.nc", &
962	& plev,play,phi,phis,presnivs, &
963	& u,v,temp,q,omega2)
964
965	CALL abort_gcm ('lmdz1d ','The End ',0)
966
967	end
968
969	#include "1DUTILS.h"
970	#include "1Dconv.h"
971
972	#endif
973

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

Download in other formats:

Original Format