Context Navigation

lmdz1d.F90 @ 2243

Last change on this file since 2243 was 2243, checked in by fhourdin, 9 years ago

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 2 :

Séparation des z0 pour la quantité de mouvement et l'enthalpie.

rugs (ou frugs, rugos, yrugos ...) disparait au profit de z0m, z0h.
Les variables qui étaient à la fois dans pbl_surface_init et

dans l'interface de pbl_surface sont suprimées de pbl_surface_init.

On travaille directement pour ces variables (evap, z0, qsol, agesno)
avec les versions de phys_state_var_mod (qui étaient
précédemment dans phys_local_var_mod

Nouveaux paramètres de contrôle :
- iflag_z0_oce (par défaut 0, et seule option active jusque là)
- z0m_seaice_omp, z0h_seaice_omp, comme leur nom l'indique (utilisées dans surf_landice
- z0min appliqué sur z0m et z0h dans pbl_surface
Introduction des fonction phyeta0_get et phyetat0_srf pour lire

les conditions de initiales dans startphy.
Du coup une seule ligne suffit pour lire et contrôler d'éventuels
problèmes.

Pour la variable fxrugs, elle est remplacée par z0m(:,nbsrf+1)

Ce choix déjà utilisé pour d'autres variables pourrait être
systématiser pour alléger l'interface de pbl_surface_mod.

Dans les sorties, les variables rugs* ont été remplacées par

des z0m* et z0h*

Nettoyage des anciens alb1/alb2 dans les lectures/écritures

des états de redémarrage (et dans pbl_surface_mod.F90).

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

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

Download in other formats:

Original Format