Context Navigation

lmdz1d.F90 @ 2019

Last change on this file since 2019 was 2019, checked in by fhourdin, 10 years ago
Passage au format libre pour inclure les ficheirs du 1D dans lmdz1d.F90
File size: 37.1 KB

Line
1	#include "../dyn3d/mod_const_mpi.F90"
2	#include "../dyn3d_common/infotrac.F90"
3	#include "../dyn3d_common/control_mod.F90"
4	#include "../dyn3d_common/disvert.F90"
5
6
7	PROGRAM lmdz1d
8
9	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
10	use phys_state_var_mod
11	use comgeomphy
12	use dimphy
13	use surface_data, only : type_ocean,ok_veget
14	use pbl_surface_mod, only : ftsoil, pbl_surface_init, &
15	& pbl_surface_final
16	use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
17
18	use infotrac ! new
19	use control_mod
20	USE indice_sol_mod
21	USE phyaqua_mod
22
23	implicit none
24	#include "dimensions.h"
25	#include "YOMCST.h"
26	#include "temps.h"
27	!!#include "control.h"
28	#include "iniprint.h"
29	#include "clesphys.h"
30	#include "dimsoil.h"
31	!#include "indicesol.h"
32
33	#include "comvert.h"
34	#include "compar1d.h"
35	#include "flux_arp.h"
36	#include "tsoilnudge.h"
37	#include "fcg_gcssold.h"
38	!!!#include "fbforcing.h"
39
40	!=====================================================================
41	! DECLARATIONS
42	!=====================================================================
43
44	!---------------------------------------------------------------------
45	! Externals
46	!---------------------------------------------------------------------
47	external fq_sat
48	real fq_sat
49
50	!---------------------------------------------------------------------
51	! Arguments d' initialisations de la physique (USER DEFINE)
52	!---------------------------------------------------------------------
53
54	integer, parameter :: ngrid=1
55	real :: zcufi = 1.
56	real :: zcvfi = 1.
57
58	!- real :: nat_surf
59	!- logical :: ok_flux_surf
60	!- real :: fsens
61	!- real :: flat
62	!- real :: tsurf
63	!- real :: rugos
64	!- real :: qsol(1:2)
65	!- real :: qsurf
66	!- real :: psurf
67	!- real :: zsurf
68	!- real :: albedo
69	!-
70	!- real :: time = 0.
71	!- real :: time_ini
72	!- real :: xlat
73	!- real :: xlon
74	!- real :: wtsurf
75	!- real :: wqsurf
76	!- real :: restart_runoff
77	!- real :: xagesno
78	!- real :: qsolinp
79	!- real :: zpicinp
80	!-
81	real :: fnday
82	real :: day, daytime
83	real :: day1
84	real :: heure
85	integer :: jour
86	integer :: mois
87	integer :: an
88
89	!---------------------------------------------------------------------
90	! Declarations related to forcing and initial profiles
91	!---------------------------------------------------------------------
92
93	integer :: kmax = llm
94	integer llm700,nq1,nq2
95	INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
96	real timestep, frac
97	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
98	real uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
99	real ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
100	real dqtdxls(nlev_max),dqtdyls(nlev_max)
101	real dqtdtls(nlev_max),thlpcar(nlev_max)
102	real qprof(nlev_max,nqmx)
103
104	! integer :: forcing_type
105	logical :: forcing_les = .false.
106	logical :: forcing_armcu = .false.
107	logical :: forcing_rico = .false.
108	logical :: forcing_radconv = .false.
109	logical :: forcing_toga = .false.
110	logical :: forcing_twpice = .false.
111	logical :: forcing_amma = .false.
112	logical :: forcing_GCM2SCM = .false.
113	logical :: forcing_GCSSold = .false.
114	logical :: forcing_sandu = .false.
115	logical :: forcing_astex = .false.
116	logical :: forcing_fire = .false.
117	integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
118	! (cf read_tsurf1d.F)
119
120	!vertical advection computation
121	! real d_t_z(llm), d_q_z(llm)
122	! real d_t_dyn_z(llm), d_q_dyn_z(llm)
123	! real zz(llm)
124	! real zfact
125
126	!flag forcings
127	logical :: nudge_wind=.true.
128	logical :: nudge_thermo=.false.
129	logical :: cptadvw=.true.
130	!=====================================================================
131	! DECLARATIONS FOR EACH CASE
132	!=====================================================================
133	!
134	#include "1D_decl_cases.h"
135	!
136	!---------------------------------------------------------------------
137	! Declarations related to vertical discretization:
138	!---------------------------------------------------------------------
139	real :: pzero=1.e5
140	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
141	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
142
143	!---------------------------------------------------------------------
144	! Declarations related to variables
145	!---------------------------------------------------------------------
146
147	real :: phi(llm)
148	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
149	real :: rlat_rad(1),rlon_rad(1)
150	real :: omega(llm+1),omega2(llm),rho(llm+1)
151	real :: ug(llm),vg(llm),fcoriolis
152	real :: sfdt, cfdt
153	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
154	real :: dt_dyn(llm)
155	real :: dt_cooling(llm),d_th_adv(llm)
156	real :: alpha
157	real :: ttt
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	! REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
164
165	!---------------------------------------------------------------------
166	! Initialization of surface variables
167	!---------------------------------------------------------------------
168	real :: run_off_lic_0(1)
169	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
170	real :: evap(1,nbsrf),frugs(1,nbsrf)
171	real :: tsoil(1,nsoilmx,nbsrf)
172	real :: agesno(1,nbsrf)
173
174	!---------------------------------------------------------------------
175	! Call to phyredem
176	!---------------------------------------------------------------------
177	logical :: ok_writedem =.true.
178
179	!---------------------------------------------------------------------
180	! Call to physiq
181	!---------------------------------------------------------------------
182	integer, parameter :: longcles=20
183	logical :: firstcall=.true.
184	logical :: lastcall=.false.
185	real :: phis = 0.0
186	real :: clesphy0(longcles) = 0.0
187	real :: dpsrf
188
189	!---------------------------------------------------------------------
190	! Initializations of boundary conditions
191	!---------------------------------------------------------------------
192	integer, parameter :: yd = 360
193	real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
194	real :: phy_alb (yd) ! Albedo land only (old value condsurf_jyg=0.3)
195	real :: phy_sst (yd) ! SST (will not be used; cf read_tsurf1d.F)
196	real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
197	real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
198	real :: phy_ice (yd) = 0.0 ! Fraction de glace
199	real :: phy_fter(yd) = 0.0 ! Fraction de terre
200	real :: phy_foce(yd) = 0.0 ! Fraction de ocean
201	real :: phy_fsic(yd) = 0.0 ! Fraction de glace
202	real :: phy_flic(yd) = 0.0 ! Fraction de glace
203
204	!---------------------------------------------------------------------
205	! Fichiers et d'autres variables
206	!---------------------------------------------------------------------
207	integer :: 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	!Al1
232	call conf_unicol
233	!Al1 moves this gcssold var from common fcg_gcssold to
234	Turb_fcg_gcssold = xTurb_fcg_gcssold
235	! --------------------------------------------------------------------
236	close(1)
237	!Al1
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	!-----------------------------------------------------------------------
331	! Choix du calendrier
332	! -------------------
333
334	! 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	!-----------------------------------------------------------------------
352	!
353	!c Date :
354	! La date est supposee donnee sous la forme [annee, numero du jour dans
355	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
356	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
357	! Le numero du jour est dans "day". L heure est traitee separement.
358	! 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	! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
366	IF(forcing_type .EQ. 61) fnday=53100./86400.
367	! 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 = int(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	if (nqtot>nqmx) STOP'Augmenter nqmx dans lmdz1d.F'
430	allocate(q(llm,nqtot)) ; q(:,:)=0.
431	allocate(dq(llm,nqtot))
432	allocate(dq_dyn(llm,nqtot))
433	allocate(d_q_adv(llm,nqtot))
434	! allocate(d_th_adv(llm))
435
436	!
437	! No ozone climatology need be read in this pre-initialization
438	! (phys_state_var_init is called again in physiq)
439	read_climoz = 0
440	!
441	call phys_state_var_init(read_climoz)
442
443	if (ngrid.ne.klon) then
444	print*,'stop in inifis'
445	print*,'Probleme de dimensions :'
446	print*,'ngrid = ',ngrid
447	print*,'klon = ',klon
448	stop
449	endif
450	!!!=====================================================================
451	!!! Feedback forcing values for Gateaux differentiation (al1)
452	!!!=====================================================================
453	!!! Surface Planck forcing bracketing call radiation
454	!! surf_Planck = 0.
455	!! surf_Conv = 0.
456	!! write(,) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
457	!!! a mettre dans le lmdz1d.def ou autre
458	!!
459	!!
460	qsol = qsolinp
461	qsurf = fq_sat(tsurf,psurf/100.)
462	rlat=xlat
463	rlon=xlon
464	day1= day_ini
465	time=daytime-day
466	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
467	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
468
469	!
470	!! mpl et jyg le 22/08/2012 :
471	!! pour que les cas a flux de surface imposes marchent
472	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
473	fsens=-wtsurfrcpdrho(1)
474	flat=-wqsurfrlvttrho(1)
475	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
476	ENDIF
477	print*,'Flux sol ',fsens,flat
478	!! ok_flux_surf=.false.
479	!! fsens=-wtsurfrcpdrho(1)
480	!! flat=-wqsurfrlvttrho(1)
481	!!!!
482
483	! Vertical discretization and pressure levels at half and mid levels:
484
485	pa = 5e4
486	!! preff= 1.01325e5
487	preff = psurf
488	IF (ok_old_disvert) THEN
489	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
490	print *,'On utilise disvert0'
491	ELSE
492	call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig, &
493	& scaleheight)
494	print *,'On utilise disvert'
495	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
496	! Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
497	ENDIF
498	sig_s=presnivs/preff
499	plev =ap+bp*psurf
500	play = 0.5*(plev(1:llm)+plev(2:llm+1))
501	!cc zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles
502
503	IF (forcing_type .eq. 59) THEN
504	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
505	write(,) '***********************'
506	do l = 1, llm
507	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
508	if (trouve_700 .and. play(l).le.70000) then
509	llm700=l
510	print *,'llm700,play=',llm700,play(l)/100.
511	trouve_700= .false.
512	endif
513	enddo
514	write(,) '***********************'
515	ENDIF
516
517	!
518	!=====================================================================
519	! EVENTUALLY, READ FORCING DATA :
520	!=====================================================================
521
522	#include "1D_read_forc_cases.h"
523
524	if (forcing_GCM2SCM) then
525	write (,) 'forcing_GCM2SCM not yet implemented'
526	stop 'in initialization'
527	endif ! forcing_GCM2SCM
528
529	print*,'mxcalc=',mxcalc
530	print*,'zlay=',zlay(mxcalc)
531	print*,'play=',play(mxcalc)
532
533	!Al1 pour SST forced, appellé depuis ocean_forced_noice
534	ts_cur = tsurf ! SST used in read_tsurf1d
535	!=====================================================================
536	! Initialisation de la physique :
537	!=====================================================================
538
539	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
540	!
541	! day_step, iphysiq lus dans gcm.def ci-dessus
542	! timestep: calcule ci-dessous from rday et day_step
543	! ngrid=1
544	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
545	! rday: defini dans suphel.F (86400.)
546	! day_ini: lu dans run.def (dayref)
547	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
548	! airefi,zcufi,zcvfi initialises au debut de ce programme
549	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
550	day_step = float(nsplit_phys)*day_step/float(iphysiq)
551	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
552	timestep =rday/day_step
553	dtime_frcg = timestep
554	!
555	zcufi=airefi
556	zcvfi=airefi
557	!
558	rlat_rad(:)=rlat(:)*rpi/180.
559	rlon_rad(:)=rlon(:)*rpi/180.
560
561	call iniphysiq(ngrid,llm,rday,day_ini,timestep, &
562	& rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
563	print*,'apres iniphysiq'
564
565	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
566	co2_ppm= 330.0
567	solaire=1370.0
568
569	! Ecriture du startphy avant le premier appel a la physique.
570	! On le met juste avant pour avoir acces a tous les champs
571	! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
572
573	if (ok_writedem) then
574
575	!--------------------------------------------------------------------------
576	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
577	! need : qsol fder snow qsurf evap rugos agesno ftsoil
578	!--------------------------------------------------------------------------
579
580	type_ocean = "force"
581	run_off_lic_0(1) = restart_runoff
582	call fonte_neige_init(run_off_lic_0)
583
584	fder=0.
585	snsrf(1,:)=0. ! couverture de neige des sous surface
586	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
587	evap=0.
588	frugs(1,:)=rugos ! couverture de neige des sous surface
589	agesno = xagesno
590	tsoil(:,:,:)=tsurf
591	!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
592	! tsoil(1,1,1)=299.18
593	! tsoil(1,2,1)=300.08
594	! tsoil(1,3,1)=301.88
595	! tsoil(1,4,1)=305.48
596	! tsoil(1,5,1)=308.00
597	! tsoil(1,6,1)=308.00
598	! tsoil(1,7,1)=308.00
599	! tsoil(1,8,1)=308.00
600	! tsoil(1,9,1)=308.00
601	! tsoil(1,10,1)=308.00
602	! tsoil(1,11,1)=308.00
603	!-----------------------------------------------------------------------
604	call pbl_surface_init(qsol, fder, snsrf, qsurfsrf, &
605	& evap, frugs, agesno, tsoil)
606
607	!------------------ prepare limit conditions for limit.nc -----------------
608	!-- Ocean force
609
610	print*,'avant phyredem'
611	pctsrf(1,:)=0.
612	if (nat_surf.eq.0.) then
613	pctsrf(1,is_oce)=1.
614	pctsrf(1,is_ter)=0.
615	else
616	pctsrf(1,is_oce)=0.
617	pctsrf(1,is_ter)=1.
618	end if
619
620	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf &
621	& ,pctsrf(1,is_oce),pctsrf(1,is_ter)
622
623	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
624	zpic = zpicinp
625	ftsol=tsurf
626	falb1 = albedo
627	falb2 = albedo
628	rugoro=rugos
629	t_ancien(1,:)=temp(:)
630	q_ancien(1,:)=q(:,1)
631	pbl_tke=1.e-8
632
633	rain_fall=0.
634	snow_fall=0.
635	solsw=0.
636	sollw=0.
637	radsol=0.
638	rnebcon=0.
639	ratqs=0.
640	clwcon=0.
641	zmea=0.
642	zstd=0.
643	zsig=0.
644	zgam=0.
645	zval=0.
646	zthe=0.
647	sig1=0.
648	w01=0.
649	u_ancien(1,:)=u(:)
650	v_ancien(1,:)=v(:)
651
652	!------------------------------------------------------------------------
653	! Make file containing restart for the physics (startphy.nc)
654	!
655	! NB: List of the variables to be written by phyredem (via put_field):
656	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
657	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
658	! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
659	! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
660	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
661	! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
662	! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,sig1,w01
663	! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
664	!------------------------------------------------------------------------
665	!Al1 =============== restart option ==========================
666	if (.not.restart) then
667	call phyredem ("startphy.nc")
668	else
669	! (desallocations)
670	print*,'callin surf final'
671	call pbl_surface_final(qsol, fder, snsrf, qsurfsrf, &
672	& evap, frugs, agesno, tsoil)
673	print*,'after surf final'
674	CALL fonte_neige_final(run_off_lic_0)
675	endif
676
677	ok_writedem=.false.
678	print*,'apres phyredem'
679
680	endif ! ok_writedem
681
682	!------------------------------------------------------------------------
683	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
684	! --------------------------------------------------
685	! NB: List of the variables to be written in limit.nc
686	! (by writelim.F, subroutine of 1DUTILS.h):
687	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
688	! phy_fter,phy_foce,phy_flic,phy_fsic)
689	!------------------------------------------------------------------------
690	do i=1,yd
691	phy_nat(i) = nat_surf
692	phy_alb(i) = albedo
693	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
694	phy_rug(i) = rugos
695	phy_fter(i) = pctsrf(1,is_ter)
696	phy_foce(i) = pctsrf(1,is_oce)
697	phy_fsic(i) = pctsrf(1,is_sic)
698	phy_flic(i) = pctsrf(1,is_lic)
699	enddo
700
701	! fabrication de limit.nc
702	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug, &
703	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
704
705
706	call phys_state_var_end
707	!Al1
708	if (restart) then
709	print*,'call to restart dyn 1d'
710	Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs, &
711	& u,v,temp,q,omega2)
712
713	print*,'fnday,annee_ref,day_ref,day_ini', &
714	& fnday,annee_ref,day_ref,day_ini
715	!** call ymds2ju(annee_ref,mois,day_ini,heure,day)
716	day = day_ini
717	day_end = day_ini + nday
718	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
719
720	! Print out the actual date of the beginning of the simulation :
721	call ju2ymds(daytime, an, mois, jour, heure)
722	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
723
724	day = int(daytime)
725	time=daytime-day
726
727	print,'*** intialised fields from restart1dyn *****'
728	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
729	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
730	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
731	! raz for safety
732	do l=1,llm
733	dq_dyn(l,1) = 0.
734	enddo
735	endif
736	!Al1 ================ end restart =================================
737	IF (ecrit_slab_oc.eq.1) then
738	open(97,file='div_slab.dat',STATUS='UNKNOWN')
739	elseif (ecrit_slab_oc.eq.0) then
740	open(97,file='div_slab.dat',STATUS='OLD')
741	endif
742	!=====================================================================
743	! START OF THE TEMPORAL LOOP :
744	!=====================================================================
745
746	do while(it.le.nint(fnday*day_step))
747
748	if (prt_level.ge.1) then
749	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
750	& it,day,time,nint(fnday*day_step)
751	print*,'PAS DE TEMPS ',timestep
752	endif
753	!Al1 demande de restartphy.nc
754	if (it.eq.nint(fnday*day_step)) lastcall=.True.
755
756	!---------------------------------------------------------------------
757	! Interpolation of forcings in time and onto model levels
758	!---------------------------------------------------------------------
759
760	#include "1D_interp_cases.h"
761
762	if (forcing_GCM2SCM) then
763	write (,) 'forcing_GCM2SCM not yet implemented'
764	stop 'in time loop'
765	endif ! forcing_GCM2SCM
766
767	!---------------------------------------------------------------------
768	! Geopotential :
769	!---------------------------------------------------------------------
770
771	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
772	do l = 1, llm-1
773	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
774	& (play(l)-play(l+1))/(play(l)+play(l+1))
775	enddo
776
777	!---------------------------------------------------------------------
778	! Listing output for debug prt_level>=1
779	!---------------------------------------------------------------------
780	if (prt_level>=1) then
781	print *,' avant physiq : -------- day time ',day,time
782	write(,) 'firstcall,lastcall,phis', &
783	& firstcall,lastcall,phis
784	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l', &
785	& 'presniv','plev','play','phi'
786	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l, &
787	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
788	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l', &
789	& 'presniv','u','v','temp','q1','q2','omega2'
790	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l, &
791	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
792	endif
793
794	!---------------------------------------------------------------------
795	! Call physiq :
796	!---------------------------------------------------------------------
797
798	call physiq(ngrid,llm, &
799	& firstcall,lastcall, &
800	& day,time,timestep, &
801	& plev,play,phi,phis,presnivs,clesphy0, &
802	& u,v,temp,q,omega2, &
803	& du_phys,dv_phys,dt_phys,dq,dpsrf, &
804	& dudyn,PVteta)
805	firstcall=.false.
806
807	!---------------------------------------------------------------------
808	! Listing output for debug prt_level>=1
809	!---------------------------------------------------------------------
810	if (prt_level>=1) then
811	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l', &
812	& 'presniv','plev','play','phi'
813	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l, &
814	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
815	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l', &
816	& 'presniv','u','v','temp','q1','q2','omega2'
817	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l, &
818	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
819	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l', &
820	& 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
821	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l, &
822	& presnivs(l),86400du_phys(l),86400dv_phys(l), &
823	& 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
824	write(,) 'dpsrf',dpsrf
825	endif
826	!---------------------------------------------------------------------
827	! Add physical tendencies :
828	!---------------------------------------------------------------------
829
830	fcoriolis=2.sin(rpixlat/180.)*romega
831	if (forcing_radconv .or. forcing_fire) then
832	fcoriolis=0.0
833	dt_cooling=0.0
834	d_th_adv=0.0
835	d_q_adv=0.0
836	endif
837
838	if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice &
839	& .or.forcing_amma) then
840	fcoriolis=0.0 ; ug=0. ; vg=0.
841	endif
842	if(forcing_rico) then
843	dt_cooling=0.
844	endif
845
846	print*, 'fcoriolis ', fcoriolis, xlat
847
848	du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
849	dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
850
851	!!!!!!!!!!!!!!!!!!!!!!!!
852	! Geostrophic wind
853	!!!!!!!!!!!!!!!!!!!!!!!!
854	sfdt = sin(0.5fcoriolistimestep)
855	cfdt = cos(0.5fcoriolistimestep)
856	!
857	du_age(1:mxcalc)= -2.sfdt/timestep &
858	& (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - &
859	& cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
860	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
861	!
862	dv_age(1:mxcalc)= -2.sfdt/timestep &
863	& (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + &
864	& sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
865	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
866	!
867	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
868	! call writefield_phy('dv_age' ,dv_age,llm)
869	! call writefield_phy('du_age' ,du_age,llm)
870	! call writefield_phy('du_phys' ,du_phys,llm)
871	! call writefield_phy('u_tend' ,u,llm)
872	! call writefield_phy('u_g' ,ug,llm)
873	!
874	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
875	!! Increment state variables
876	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
877	! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
878	! au dessus de 700hpa, on relaxe vers les profils initiaux
879	if (forcing_sandu .OR. forcing_astex) then
880	#include "1D_nudge_sandu_astex.h"
881	else
882	u(1:mxcalc)=u(1:mxcalc) + timestep*( &
883	& du_phys(1:mxcalc) &
884	& +du_age(1:mxcalc) )
885	v(1:mxcalc)=v(1:mxcalc) + timestep*( &
886	& dv_phys(1:mxcalc) &
887	& +dv_age(1:mxcalc) )
888	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( &
889	& dq(1:mxcalc,:) &
890	& +d_q_adv(1:mxcalc,:) )
891
892	if (prt_level.ge.1) then
893	print *, &
894	& 'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ', &
895	& temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
896	print*,du_phys
897	print*, v
898	print*, vg
899	endif
900
901	temp(1:mxcalc)=temp(1:mxcalc)+timestep*( &
902	& dt_phys(1:mxcalc) &
903	& +d_th_adv(1:mxcalc) &
904	& +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
905
906	endif ! forcing_sandu or forcing_astex
907
908	teta=temp(pzero/play)*rkappa
909	!
910	!---------------------------------------------------------------------
911	! Nudge soil temperature if requested
912	!---------------------------------------------------------------------
913
914	IF (nudge_tsoil) THEN
915	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) &
916	& -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
917	ENDIF
918
919	!---------------------------------------------------------------------
920	! Add large-scale tendencies (advection, etc) :
921	!---------------------------------------------------------------------
922
923	!cc nrlmd
924	!cc tmpvar=teta
925	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
926	!cc
927	!cc teta(1:mxcalc)=tmpvar(1:mxcalc)
928	!cc tmpvar(:)=q(:,1)
929	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
930	!cc q(1:mxcalc,1)=tmpvar(1:mxcalc)
931	!cc tmpvar(:)=q(:,2)
932	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
933	!cc q(1:mxcalc,2)=tmpvar(1:mxcalc)
934
935	!---------------------------------------------------------------------
936	! Air temperature :
937	!---------------------------------------------------------------------
938	if (lastcall) then
939	print*,'Pas de temps final ',it
940	call ju2ymds(daytime, an, mois, jour, heure)
941	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
942	endif
943
944	! incremente day time
945	! print*,'daytime bef',daytime,1./day_step
946	daytime = daytime+1./day_step
947	!Al1dbg
948	day = int(daytime+0.1/day_step)
949	! time = max(daytime-day,0.0)
950	!Al1&jyg: correction de bug
951	!cc time = real(mod(it,day_step))/day_step
952	time = time_ini/24.+real(mod(it,day_step))/day_step
953	! print*,'daytime nxt time',daytime,time
954	it=it+1
955
956	enddo
957
958	!Al1
959	if (ecrit_slab_oc.ne.-1) close(97)
960
961	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
962	! -------------------------------------
963	call dyn1dredem("restart1dyn.nc", &
964	& plev,play,phi,phis,presnivs, &
965	& u,v,temp,q,omega2)
966
967	CALL abort_gcm ('lmdz1d ','The End ',0)
968
969	end
970
971	#include "1DUTILS.h"
972	#include "1Dconv.h"
973
974

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

Download in other formats:

Original Format