Context Navigation

lmdz1d.F @ 1960

Last change on this file since 1960 was 1960, checked in by fhourdin, 11 years ago

Nettoyage du code 1D pour limiter les warning en mode debug.
Essentiellement des declarations inutiles.

Cleaning of the 1D code in order to limit warnings in debug mode.
Essentially supressing unused declared variables.

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

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

Download in other formats:

Original Format