Context Navigation

scm.F90 @ 3890

Last change on this file since 3890 was 3888, checked in by jyg, 4 years ago

New provisional version of the splitting of the
diffusive boundary layer into inwake and offwake
PBLs. The splitting of the diffuse BL should NOT
be activated yet for general purpose simulations.

The splitting is activated by:
mod(iflag_pbl_split,10)=1 for the option with
fixed surface temperature and
mod(iflag_pbl_split,10)=2 for the option with
coupled surface temperature.

iflag_pbl_split=0 ==> no splittingat all.
iflag_pbl_split=10 ==> splitting of thermals.
iflag_pbl_split=11 ==> splitting of thermals and
of vertical diffusion (fixed surf. temp.).
iflag_pbl_split=12 ==> splitting of thermals and
of vertical diffusion (coupled surf. temp.).

File size: 39.5 KB

Line
1	SUBROUTINE scm
2
3	USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
4	USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
5	clwcon, detr_therm, &
6	qsol, fevap, z0m, z0h, agesno, &
7	du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
8	falb_dir, falb_dif, &
9	ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
10	rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
11	solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
12	wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
13	wake_deltaq, wake_deltat, wake_s, wake_dens, &
14	zgam, zmax0, zmea, zpic, zsig, &
15	zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, ql_ancien, qs_ancien, &
16	prlw_ancien, prsw_ancien, prw_ancien, &
17	u10m,v10m,ale_wake,ale_bl_stat
18
19
20	USE dimphy
21	USE surface_data, only : type_ocean,ok_veget
22	USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
23	pbl_surface_final
24	USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
25
26	USE infotrac ! new
27	USE control_mod
28	USE indice_sol_mod
29	USE phyaqua_mod
30	! USE mod_1D_cases_read
31	USE mod_1D_cases_read_std
32	!USE mod_1D_amma_read
33	USE print_control_mod, ONLY: lunout, prt_level
34	USE iniphysiq_mod, ONLY: iniphysiq
35	USE mod_const_mpi, ONLY: comm_lmdz
36	USE physiq_mod, ONLY: physiq
37	USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
38	preff, aps, bps, pseudoalt, scaleheight
39	USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
40	itau_dyn, itau_phy, start_time, year_len
41	USE phys_cal_mod, ONLY : year_len_phys_cal_mod => year_len
42
43	implicit none
44	#include "dimensions.h"
45	#include "YOMCST.h"
46	!!#include "control.h"
47	#include "clesphys.h"
48	#include "dimsoil.h"
49	!#include "indicesol.h"
50
51	#include "compar1d.h"
52	#include "flux_arp.h"
53	#include "date_cas.h"
54	#include "tsoilnudge.h"
55	#include "fcg_gcssold.h"
56	#include "compbl.h"
57
58	!=====================================================================
59	! DECLARATIONS
60	!=====================================================================
61
62	#undef OUTPUT_PHYS_SCM
63
64	!---------------------------------------------------------------------
65	! Externals
66	!---------------------------------------------------------------------
67	external fq_sat
68	real fq_sat
69
70	!---------------------------------------------------------------------
71	! Arguments d' initialisations de la physique (USER DEFINE)
72	!---------------------------------------------------------------------
73
74	integer, parameter :: ngrid=1
75	real :: zcufi = 1.
76	real :: zcvfi = 1.
77	real :: fnday
78	real :: day, daytime
79	real :: day1
80	real :: heure
81	integer :: jour
82	integer :: mois
83	integer :: an
84
85	!---------------------------------------------------------------------
86	! Declarations related to forcing and initial profiles
87	!---------------------------------------------------------------------
88
89	integer :: kmax = llm
90	integer llm700,nq1,nq2
91	INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
92	real timestep, frac
93	real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
94	real uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
95	real ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
96	real dqtdxls(nlev_max),dqtdyls(nlev_max)
97	real dqtdtls(nlev_max),thlpcar(nlev_max)
98	real qprof(nlev_max,nqmx)
99
100	! integer :: forcing_type
101	logical :: forcing_les = .false.
102	logical :: forcing_armcu = .false.
103	logical :: forcing_rico = .false.
104	logical :: forcing_radconv = .false.
105	logical :: forcing_toga = .false.
106	logical :: forcing_twpice = .false.
107	logical :: forcing_amma = .false.
108	logical :: forcing_dice = .false.
109	logical :: forcing_gabls4 = .false.
110
111	logical :: forcing_GCM2SCM = .false.
112	logical :: forcing_GCSSold = .false.
113	logical :: forcing_sandu = .false.
114	logical :: forcing_astex = .false.
115	logical :: forcing_fire = .false.
116	logical :: forcing_case = .false.
117	logical :: forcing_case2 = .false.
118	logical :: forcing_SCM = .false.
119
120	!flag forcings
121	logical :: nudge_wind=.true.
122	logical :: nudge_thermo=.false.
123	logical :: cptadvw=.true.
124
125
126	!=====================================================================
127	! DECLARATIONS FOR EACH CASE
128	!=====================================================================
129	!
130	#include "1D_decl_cases.h"
131	!
132	!---------------------------------------------------------------------
133	! Declarations related to nudging
134	!---------------------------------------------------------------------
135	integer :: nudge_max
136	parameter (nudge_max=9)
137	integer :: inudge_RHT=1
138	integer :: inudge_UV=2
139	logical :: nudge(nudge_max)
140	real :: t_targ(llm)
141	real :: rh_targ(llm)
142	real :: u_targ(llm)
143	real :: v_targ(llm)
144	!
145	!---------------------------------------------------------------------
146	! Declarations related to vertical discretization:
147	!---------------------------------------------------------------------
148	real :: pzero=1.e5
149	real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
150	real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
151
152	!---------------------------------------------------------------------
153	! Declarations related to variables
154	!---------------------------------------------------------------------
155
156	real :: phi(llm)
157	real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
158	REAL rot(1, llm) ! relative vorticity, in s-1
159	real :: rlat_rad(1),rlon_rad(1)
160	real :: omega(llm),omega2(llm),rho(llm+1)
161	real :: ug(llm),vg(llm),fcoriolis
162	real :: sfdt, cfdt
163	real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
164	real :: w_adv(llm),z_adv(llm)
165	real :: d_t_vert_adv(llm),d_u_vert_adv(llm),d_v_vert_adv(llm)
166	real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
167	real :: d_u_nudge(llm),d_v_nudge(llm)
168	! real :: d_u_adv(llm),d_v_adv(llm)
169	real :: d_u_age(llm),d_v_age(llm)
170	real :: alpha
171	real :: ttt
172
173	REAL, ALLOCATABLE, DIMENSION(:,:):: q
174	REAL, ALLOCATABLE, DIMENSION(:,:):: dq
175	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_vert_adv
176	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
177	REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
178	! REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
179
180	!---------------------------------------------------------------------
181	! Initialization of surface variables
182	!---------------------------------------------------------------------
183	real :: run_off_lic_0(1)
184	real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
185	real :: tsoil(1,nsoilmx,nbsrf)
186	! real :: agesno(1,nbsrf)
187
188	!---------------------------------------------------------------------
189	! Call to phyredem
190	!---------------------------------------------------------------------
191	logical :: ok_writedem =.true.
192	real :: sollw_in = 0.
193	real :: solsw_in = 0.
194
195	!---------------------------------------------------------------------
196	! Call to physiq
197	!---------------------------------------------------------------------
198	logical :: firstcall=.true.
199	logical :: lastcall=.false.
200	real :: phis(1) = 0.0
201	real :: dpsrf(1)
202
203	!---------------------------------------------------------------------
204	! Initializations of boundary conditions
205	!---------------------------------------------------------------------
206	real, allocatable :: phy_nat (:) ! 0=ocean libre,1=land,2=glacier,3=banquise
207	real, allocatable :: phy_alb (:) ! Albedo land only (old value condsurf_jyg=0.3)
208	real, allocatable :: phy_sst (:) ! SST (will not be used; cf read_tsurf1d.F)
209	real, allocatable :: phy_bil (:) ! Ne sert que pour les slab_ocean
210	real, allocatable :: phy_rug (:) ! Longueur rugosite utilisee sur land only
211	real, allocatable :: phy_ice (:) ! Fraction de glace
212	real, allocatable :: phy_fter(:) ! Fraction de terre
213	real, allocatable :: phy_foce(:) ! Fraction de ocean
214	real, allocatable :: phy_fsic(:) ! Fraction de glace
215	real, allocatable :: phy_flic(:) ! Fraction de glace
216
217	!---------------------------------------------------------------------
218	! Fichiers et d'autres variables
219	!---------------------------------------------------------------------
220	integer :: k,l,i,it=1,mxcalc
221	integer :: nsrf
222	integer jcode
223	INTEGER read_climoz
224	!
225	integer :: it_end ! iteration number of the last call
226	!Al1,plev,play,phi,phis,presnivs,
227	integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
228	data ecrit_slab_oc/-1/
229	!
230	! if flag_inhib_forcing = 0, tendencies of forcing are added
231	! <> 0, tendencies of forcing are not added
232	INTEGER :: flag_inhib_forcing = 0
233
234
235	print*,'VOUS ENTREZ DANS LE 1D FORMAT STANDARD'
236
237	!=====================================================================
238	! INITIALIZATIONS
239	!=====================================================================
240	du_phys(:)=0.
241	dv_phys(:)=0.
242	dt_phys(:)=0.
243	d_t_vert_adv(:)=0.
244	d_u_vert_adv(:)=0.
245	d_v_vert_adv(:)=0.
246	dt_cooling(:)=0.
247	d_t_adv(:)=0.
248	d_t_nudge(:)=0.
249	d_u_nudge(:)=0.
250	d_v_nudge(:)=0.
251	d_u_adv(:)=0.
252	d_v_adv(:)=0.
253	d_u_age(:)=0.
254	d_v_age(:)=0.
255
256
257	! Initialization of Common turb_forcing
258	dtime_frcg = 0.
259	Turb_fcg_gcssold=.false.
260	hthturb_gcssold = 0.
261	hqturb_gcssold = 0.
262
263
264
265
266	!---------------------------------------------------------------------
267	! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
268	!---------------------------------------------------------------------
269	call conf_unicol
270	!Al1 moves this gcssold var from common fcg_gcssold to
271	Turb_fcg_gcssold = xTurb_fcg_gcssold
272	! --------------------------------------------------------------------
273	close(1)
274	write(,) 'lmdz1d.def lu => unicol.def'
275
276	forcing_SCM = .true.
277	year_ini_cas=1997
278	! It is possible that those parameters are run twice.
279	! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
280
281
282	call getin('anneeref',year_ini_cas)
283	call getin('dayref',day_deb)
284	mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
285	call getin('time_ini',heure_ini_cas)
286
287	print*,'NATURE DE LA SURFACE ',nat_surf
288	!
289	! Initialization of the logical switch for nudging
290
291	jcode = iflag_nudge
292	do i = 1,nudge_max
293	nudge(i) = mod(jcode,10) .ge. 1
294	jcode = jcode/10
295	enddo
296	!-----------------------------------------------------------------------
297	! Definition of the run
298	!-----------------------------------------------------------------------
299
300	call conf_gcm( 99, .TRUE. )
301
302	!-----------------------------------------------------------------------
303	allocate( phy_nat (year_len)) ! 0=ocean libre,1=land,2=glacier,3=banquise
304	phy_nat(:)=0.0
305	allocate( phy_alb (year_len)) ! Albedo land only (old value condsurf_jyg=0.3)
306	allocate( phy_sst (year_len)) ! SST (will not be used; cf read_tsurf1d.F)
307	allocate( phy_bil (year_len)) ! Ne sert que pour les slab_ocean
308	phy_bil(:)=1.0
309	allocate( phy_rug (year_len)) ! Longueur rugosite utilisee sur land only
310	allocate( phy_ice (year_len)) ! Fraction de glace
311	phy_ice(:)=0.0
312	allocate( phy_fter(year_len)) ! Fraction de terre
313	phy_fter(:)=0.0
314	allocate( phy_foce(year_len)) ! Fraction de ocean
315	phy_foce(:)=0.0
316	allocate( phy_fsic(year_len)) ! Fraction de glace
317	phy_fsic(:)=0.0
318	allocate( phy_flic(year_len)) ! Fraction de glace
319	phy_flic(:)=0.0
320
321
322	!-----------------------------------------------------------------------
323	! Choix du calendrier
324	! -------------------
325
326	! 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	stop 'gregorian calend should not be used by normal user'
338	call ioconf_calendar('gregorian') ! not to be used by normal users
339	write(,)'CALENDRIER CHOISI: Gregorien'
340	else
341	write (,) 'ERROR : unknown calendar ', calend
342	stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
343	endif
344	!-----------------------------------------------------------------------
345	!
346	!c Date :
347	! La date est supposee donnee sous la forme [annee, numero du jour dans
348	! l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
349	! On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
350	! Le numero du jour est dans "day". L heure est traitee separement.
351	! La date complete est dans "daytime" (l'unite est le jour).
352
353
354	if (nday>0) then
355	fnday=nday
356	else
357	fnday=-nday/float(day_step)
358	endif
359	print *,'fnday=',fnday
360	! start_time doit etre en FRACTION DE JOUR
361	start_time=time_ini/24.
362
363	annee_ref = anneeref
364	mois = 1
365	day_ref = dayref
366	heure = 0.
367	itau_dyn = 0
368	itau_phy = 0
369	call ymds2ju(annee_ref,mois,day_ref,heure,day)
370	day_ini = int(day)
371	day_end = day_ini + int(fnday)
372
373	! Convert the initial date to Julian day
374	day_ini_cas=day_deb
375	print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas
376	call ymds2ju &
377	& (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600 &
378	& ,day_ju_ini_cas)
379	print*,'time case 2',day_ini_cas,day_ju_ini_cas
380	daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation
381
382	! Print out the actual date of the beginning of the simulation :
383	call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
384	print *,' Time of beginning : ', &
385	& year_print, month_print, day_print, sec_print
386
387	!---------------------------------------------------------------------
388	! Initialization of dimensions, geometry and initial state
389	!---------------------------------------------------------------------
390	! call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
391	! but we still need to initialize dimphy module (klon,klev,etc.) here.
392	call init_dimphy1D(1,llm)
393	call suphel
394	call infotrac_init
395
396	if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
397	allocate(q(llm,nqtot)) ; q(:,:)=0.
398	allocate(dq(llm,nqtot))
399	allocate(d_q_vert_adv(llm,nqtot))
400	allocate(d_q_adv(llm,nqtot))
401	allocate(d_q_nudge(llm,nqtot))
402	! allocate(d_th_adv(llm))
403
404	q(:,:) = 0.
405	dq(:,:) = 0.
406	d_q_vert_adv(:,:) = 0.
407	d_q_adv(:,:) = 0.
408	d_q_nudge(:,:) = 0.
409
410	!
411	! No ozone climatology need be read in this pre-initialization
412	! (phys_state_var_init is called again in physiq)
413	read_climoz = 0
414	nsw=6
415
416	call phys_state_var_init(read_climoz)
417
418	if (ngrid.ne.klon) then
419	print*,'stop in inifis'
420	print*,'Probleme de dimensions :'
421	print*,'ngrid = ',ngrid
422	print*,'klon = ',klon
423	stop
424	endif
425	!!!=====================================================================
426	!!! Feedback forcing values for Gateaux differentiation (al1)
427	!!!=====================================================================
428	!!
429	qsol = qsolinp
430	qsurf = fq_sat(tsurf,psurf/100.)
431	beta_aridity(:,:) = beta_surf
432	day1= day_ini
433	time=daytime-day
434	ts_toga(1)=tsurf ! needed by read_tsurf1d.F
435	rho(1)=psurf/(rdtsurf(1.+(rv/rd-1.)*qsurf))
436
437	!
438	!! mpl et jyg le 22/08/2012 :
439	!! pour que les cas a flux de surface imposes marchent
440	IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
441	fsens=-wtsurfrcpdrho(1)
442	flat=-wqsurfrlvttrho(1)
443	print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
444	ENDIF
445	print*,'Flux sol ',fsens,flat
446
447	! Vertical discretization and pressure levels at half and mid levels:
448
449	pa = 5e4
450	!! preff= 1.01325e5
451	preff = psurf
452	IF (ok_old_disvert) THEN
453	call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
454	print *,'On utilise disvert0'
455	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
456	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
457	scaleheight=8.
458	pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
459	ELSE
460	call disvert()
461	print *,'On utilise disvert'
462	! Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
463	! Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
464	ENDIF
465
466	sig_s=presnivs/preff
467	plev =ap+bp*psurf
468	play = 0.5*(plev(1:llm)+plev(2:llm+1))
469	zlay=-rd300.log(play/psurf)/rg ! moved after reading profiles.
470
471	IF (forcing_type .eq. 59) THEN
472	! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
473	write(,) '***********************'
474	do l = 1, llm
475	write(,) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
476	if (trouve_700 .and. play(l).le.70000) then
477	llm700=l
478	print *,'llm700,play=',llm700,play(l)/100.
479	trouve_700= .false.
480	endif
481	enddo
482	write(,) '***********************'
483	ENDIF
484
485	!
486	!=====================================================================
487	! EVENTUALLY, READ FORCING DATA :
488	!=====================================================================
489
490	#include "1D_read_forc_cases.h"
491	print,'A d_t_adv ',d_t_adv(1:20)86400
492
493	if (forcing_GCM2SCM) then
494	write (,) 'forcing_GCM2SCM not yet implemented'
495	stop 'in initialization'
496	endif ! forcing_GCM2SCM
497
498	print*,'mxcalc=',mxcalc
499	! print*,'zlay=',zlay(mxcalc)
500	! print*,'play=',play(mxcalc)
501
502	!! When surface temperature is forced
503	tg= tsurf ! surface T used in read_tsurf1d
504
505
506	!=====================================================================
507	! Initialisation de la physique :
508	!=====================================================================
509
510	! Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
511	!
512	! day_step, iphysiq lus dans gcm.def ci-dessus
513	! timestep: calcule ci-dessous from rday et day_step
514	! ngrid=1
515	! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
516	! rday: defini dans suphel.F (86400.)
517	! day_ini: lu dans run.def (dayref)
518	! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
519	! airefi,zcufi,zcvfi initialises au debut de ce programme
520	! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
521
522
523	day_step = float(nsplit_phys)*day_step/float(iphysiq)
524	write (,) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
525	timestep =rday/day_step
526	dtime_frcg = timestep
527	!
528	zcufi=airefi
529	zcvfi=airefi
530	!
531	rlat_rad(1)=xlat*rpi/180.
532	rlon_rad(1)=xlon*rpi/180.
533
534	! iniphysiq will call iniaqua who needs year_len from phys_cal_mod
535	year_len_phys_cal_mod=year_len
536
537	! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
538	! e.g. for cell boundaries, which are meaningless in 1D; so pad these
539	! with '0.' when necessary
540
541	call iniphysiq(iim,jjm,llm, &
542	1,comm_lmdz, &
543	rday,day_ini,timestep, &
544	(/rlat_rad(1),0./),(/0./), &
545	(/0.,0./),(/rlon_rad(1),0./), &
546	(/ (/airefi,0./),(/0.,0./) /), &
547	(/zcufi,0.,0.,0./), &
548	(/zcvfi,0./), &
549	ra,rg,rd,rcpd,1)
550	print*,'apres iniphysiq'
551
552	! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
553	co2_ppm= 330.0
554	solaire=1370.0
555
556	! Ecriture du startphy avant le premier appel a la physique.
557	! On le met juste avant pour avoir acces a tous les champs
558
559	if (ok_writedem) then
560
561	!--------------------------------------------------------------------------
562	! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
563	! need : qsol fder snow qsurf evap rugos agesno ftsoil
564	!--------------------------------------------------------------------------
565
566	type_ocean = "force"
567	run_off_lic_0(1) = restart_runoff
568	call fonte_neige_init(run_off_lic_0)
569
570	fder=0.
571	snsrf(1,:)=snowmass ! masse de neige des sous surface
572	qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
573	fevap=0.
574	z0m(1,:)=rugos ! couverture de neige des sous surface
575	z0h(1,:)=rugosh ! couverture de neige des sous surface
576	agesno = xagesno
577	tsoil(:,:,:)=tsurf
578	!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
579	! tsoil(1,1,1)=299.18
580	! tsoil(1,2,1)=300.08
581	! tsoil(1,3,1)=301.88
582	! tsoil(1,4,1)=305.48
583	! tsoil(1,5,1)=308.00
584	! tsoil(1,6,1)=308.00
585	! tsoil(1,7,1)=308.00
586	! tsoil(1,8,1)=308.00
587	! tsoil(1,9,1)=308.00
588	! tsoil(1,10,1)=308.00
589	! tsoil(1,11,1)=308.00
590	!-----------------------------------------------------------------------
591	call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
592
593	!------------------ prepare limit conditions for limit.nc -----------------
594	!-- Ocean force
595
596	print*,'avant phyredem'
597	pctsrf(1,:)=0.
598	if (nat_surf.eq.0.) then
599	pctsrf(1,is_oce)=1.
600	pctsrf(1,is_ter)=0.
601	pctsrf(1,is_lic)=0.
602	pctsrf(1,is_sic)=0.
603	else if (nat_surf .eq. 1) then
604	pctsrf(1,is_oce)=0.
605	pctsrf(1,is_ter)=1.
606	pctsrf(1,is_lic)=0.
607	pctsrf(1,is_sic)=0.
608	else if (nat_surf .eq. 2) then
609	pctsrf(1,is_oce)=0.
610	pctsrf(1,is_ter)=0.
611	pctsrf(1,is_lic)=1.
612	pctsrf(1,is_sic)=0.
613	else if (nat_surf .eq. 3) then
614	pctsrf(1,is_oce)=0.
615	pctsrf(1,is_ter)=0.
616	pctsrf(1,is_lic)=0.
617	pctsrf(1,is_sic)=1.
618
619	end if
620
621
622	print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf &
623	& ,pctsrf(1,is_oce),pctsrf(1,is_ter)
624
625	zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
626	zpic = zpicinp
627	ftsol=tsurf
628	falb_dir=albedo
629	falb_dif=albedo
630	rugoro=rugos
631	t_ancien(1,:)=temp(:)
632	q_ancien(1,:)=q(:,1)
633	ql_ancien = 0.
634	qs_ancien = 0.
635	prlw_ancien = 0.
636	prsw_ancien = 0.
637	prw_ancien = 0.
638	!jyg<
639	! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
640	!! pbl_tke(:,:,:)=1.e-8
641	! pbl_tke(:,:,:)=0.
642	! pbl_tke(:,2,:)=1.e-2
643	!>jyg
644	rain_fall=0.
645	snow_fall=0.
646	solsw=0.
647	sollw=0.
648	sollwdown=rsigmatsurf*4
649	radsol=0.
650	rnebcon=0.
651	ratqs=0.
652	clwcon=0.
653	zmax0 = 0.
654	zmea=0.
655	zstd=0.
656	zsig=0.
657	zgam=0.
658	zval=0.
659	zthe=0.
660	sig1=0.
661	w01=0.
662	wake_cstar = 0.
663	wake_deltaq = 0.
664	wake_deltat = 0.
665	wake_delta_pbl_TKE(:,:,:) = 0.
666	delta_tsurf = 0.
667	wake_fip = 0.
668	wake_pe = 0.
669	wake_s = 0.
670	wake_dens = 0.
671	ale_bl = 0.
672	ale_bl_trig = 0.
673	alp_bl = 0.
674	IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
675	IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
676	entr_therm = 0.
677	detr_therm = 0.
678	f0 = 0.
679	fm_therm = 0.
680	u_ancien(1,:)=u(:)
681	v_ancien(1,:)=v(:)
682
683	u10m=0.
684	v10m=0.
685	ale_wake=0.
686	ale_bl_stat=0.
687
688	!------------------------------------------------------------------------
689	! Make file containing restart for the physics (startphy.nc)
690	!
691	! NB: List of the variables to be written by phyredem (via put_field):
692	! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
693	! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
694	! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
695	! radsol,solsw,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
696	! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
697	! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
698	! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
699	! wake_deltat,wake_deltaq,wake_s,wake_dens,wake_cstar,
700	! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
701	!
702	! NB2: The content of the startphy.nc file depends on some flags defined in
703	! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
704	! to be set at some arbitratry convenient values.
705	!------------------------------------------------------------------------
706	!Al1 =============== restart option ======================================
707	if (.not.restart) then
708	iflag_pbl = 5
709	call phyredem ("startphy.nc")
710	else
711	! (desallocations)
712	print*,'callin surf final'
713	call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
714	print*,'after surf final'
715	CALL fonte_neige_final(run_off_lic_0)
716	endif
717
718	ok_writedem=.false.
719	print*,'apres phyredem'
720
721	endif ! ok_writedem
722
723	!------------------------------------------------------------------------
724	! Make file containing boundary conditions (limit.nc) Al1->restartdyn*
725	! --------------------------------------------------
726	! NB: List of the variables to be written in limit.nc
727	! (by writelim.F, subroutine of 1DUTILS.h):
728	! phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
729	! phy_fter,phy_foce,phy_flic,phy_fsic)
730	!------------------------------------------------------------------------
731	do i=1,year_len
732	phy_nat(i) = nat_surf
733	phy_alb(i) = albedo
734	phy_sst(i) = tsurf ! read_tsurf1d will be used instead
735	phy_rug(i) = rugos
736	phy_fter(i) = pctsrf(1,is_ter)
737	phy_foce(i) = pctsrf(1,is_oce)
738	phy_fsic(i) = pctsrf(1,is_sic)
739	phy_flic(i) = pctsrf(1,is_lic)
740	enddo
741
742	! fabrication de limit.nc
743	call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug, &
744	& phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
745
746
747	call phys_state_var_end
748	!Al1
749	if (restart) then
750	print*,'call to restart dyn 1d'
751	Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs, &
752	& u,v,temp,q,omega2)
753
754	print*,'fnday,annee_ref,day_ref,day_ini', &
755	& fnday,annee_ref,day_ref,day_ini
756	!** call ymds2ju(annee_ref,mois,day_ini,heure,day)
757	day = day_ini
758	day_end = day_ini + nday
759	daytime = day + time_ini/24. ! 1st day and initial time of the simulation
760
761	! Print out the actual date of the beginning of the simulation :
762	call ju2ymds(daytime, an, mois, jour, heure)
763	print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
764
765	day = int(daytime)
766	time=daytime-day
767
768	print,'*** intialised fields from restart1dyn *****'
769	print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
770	print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
771	print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
772	! raz for safety
773	do l=1,llm
774	d_q_vert_adv(l,1) = 0.
775	enddo
776	endif
777	!====================== end restart =================================
778	IF (ecrit_slab_oc.eq.1) then
779	open(97,file='div_slab.dat',STATUS='UNKNOWN')
780	elseif (ecrit_slab_oc.eq.0) then
781	open(97,file='div_slab.dat',STATUS='OLD')
782	endif
783	!
784	!=====================================================================
785	#ifdef OUTPUT_PHYS_SCM
786	CALL iophys_ini
787	#endif
788
789	!=====================================================================
790	! START OF THE TEMPORAL LOOP :
791	!=====================================================================
792
793	it_end = nint(fnday*day_step)
794	do while(it.le.it_end)
795
796	if (prt_level.ge.1) then
797	print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', &
798	& it,day,time,it_end,day_step
799	print*,'PAS DE TEMPS ',timestep
800	endif
801	if (it.eq.it_end) lastcall=.True.
802
803	!---------------------------------------------------------------------
804	! Interpolation of forcings in time and onto model levels
805	!---------------------------------------------------------------------
806
807	#include "1D_interp_cases.h"
808
809	!---------------------------------------------------------------------
810	! Geopotential :
811	!---------------------------------------------------------------------
812	! phis(1)=zsurf*RG
813	! phi(1)=phis(1)+RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
814	phi(1)=RDtemp(1)(plev(1)-play(1))/(.5*(plev(1)+play(1)))
815
816	do l = 1, llm-1
817	phi(l+1)=phi(l)+RD(temp(l)+temp(l+1)) &
818	& (play(l)-play(l+1))/(play(l)+play(l+1))
819	enddo
820
821
822	!---------------------------------------------------------------------
823	! Vertical advection
824	!---------------------------------------------------------------------
825
826	IF ( forc_w+forc_omega > 0 ) THEN
827
828	IF ( forc_w == 1 ) THEN
829	w_adv=w_mod_cas
830	z_adv=phi/RG
831	ELSE
832	w_adv=omega
833	z_adv=play
834	ENDIF
835
836	teta=temp(pzero/play)*rkappa
837	do l=2,llm-1
838	! vertical tendencies computed as d X / d t = -W d X / d z
839	d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
840	d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
841	! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
842	d_t_vert_adv(l)=-w_adv(l)(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))*rkappa
843	d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))
844	enddo
845	d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
846	d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
847	d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
848	d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
849
850	print*,'OMEGA ',w_adv(10),z_adv(10)
851
852	ENDIF
853
854	!---------------------------------------------------------------------
855	! Listing output for debug prt_level>=1
856	!---------------------------------------------------------------------
857	if (prt_level>=1) then
858	print *,' avant physiq : -------- day time ',day,time
859	write(,) 'firstcall,lastcall,phis', &
860	& firstcall,lastcall,phis
861	end if
862	if (prt_level>=5) then
863	write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l', &
864	& 'presniv','plev','play','phi'
865	write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l, &
866	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
867	write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l', &
868	& 'presniv','u','v','temp','q1','q2','omega2'
869	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l, &
870	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
871	endif
872
873	!---------------------------------------------------------------------
874	! Call physiq :
875	!---------------------------------------------------------------------
876	call physiq(ngrid,llm, &
877	firstcall,lastcall,timestep, &
878	plev,play,phi,phis,presnivs, &
879	u,v, rot, temp,q,omega2, &
880	du_phys,dv_phys,dt_phys,dq,dpsrf)
881	firstcall=.false.
882
883	!---------------------------------------------------------------------
884	! Listing output for debug
885	!---------------------------------------------------------------------
886	if (prt_level>=5) then
887	write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l', &
888	& 'presniv','plev','play','phi'
889	write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l, &
890	& presnivs(l),plev(l),play(l),phi(l),l=1,llm)
891	write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l', &
892	& 'presniv','u','v','temp','q1','q2','omega2'
893	write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l, &
894	& presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
895	write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l', &
896	& 'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'
897	write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l, &
898	& presnivs(l),86400du_phys(l),86400dv_phys(l), &
899	& 86400dt_phys(l),86400dq(l,1),dq(l,2),l=1,llm)
900	write(,) 'dpsrf',dpsrf
901	endif
902	!---------------------------------------------------------------------
903	! Add physical tendencies :
904	!---------------------------------------------------------------------
905
906	fcoriolis=2.sin(rpixlat/180.)*romega
907
908	IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
909	fcoriolis, xlat,mxcalc
910
911	!---------------------------------------------------------------------
912	! Geostrophic forcing
913	!---------------------------------------------------------------------
914
915	IF ( forc_geo == 0 ) THEN
916	d_u_age(1:mxcalc)=0.
917	d_v_age(1:mxcalc)=0.
918	ELSE
919	sfdt = sin(0.5fcoriolistimestep)
920	cfdt = cos(0.5fcoriolistimestep)
921
922	d_u_age(1:mxcalc)= -2.sfdt/timestep &
923	& (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - &
924	& cfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
925	!! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
926	!
927	d_v_age(1:mxcalc)= -2.sfdt/timestep &
928	& (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + &
929	& sfdt*(v(1:mxcalc)-vg(1:mxcalc)) )
930	!! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
931	ENDIF
932	!
933	!---------------------------------------------------------------------
934	! Nudging
935	!---------------------------------------------------------------------
936	d_t_nudge(:) = 0.
937	d_u_nudge(:) = 0.
938	d_v_nudge(:) = 0.
939	d_q_nudge(:,:) = 0.
940	DO l=1,llm
941	IF ( play(l) < p_nudging_u .AND. nint(nudging_u) /= 0 ) &
942	& d_u_nudge(l)=(u_nudg_mod_cas(l)-u(l))/nudging_u
943	IF ( play(l) < p_nudging_v .AND. nint(nudging_v) /= 0 ) &
944	& d_v_nudge(l)=(v_nudg_mod_cas(l)-v(l))/nudging_v
945	IF ( play(l) < p_nudging_t .AND. nint(nudging_t) /= 0 ) &
946	& d_t_nudge(l)=(temp_nudg_mod_cas(l)-temp(l))/nudging_t
947	IF ( play(l) < p_nudging_qv .AND. nint(nudging_qv) /= 0 ) &
948	& d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv
949	ENDDO
950
951	!---------------------------------------------------------------------
952	! Optional outputs
953	!---------------------------------------------------------------------
954	#ifdef OUTPUT_PHYS_SCM
955	CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
956	CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv)
957	CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv)
958	CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv)
959	CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1))
960	CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas)
961	CALL iophys_ecrit('u',klev,'u','m/s',u)
962	CALL iophys_ecrit('unud',klev,'unud','m/s',u_nudg_mod_cas)
963	CALL iophys_ecrit('v',klev,'v','m/s',v)
964	CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas)
965	CALL iophys_ecrit('temp',klev,'temp','K',temp)
966	CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas)
967	CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge)
968	CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',100086400d_q_nudge(:,1))
969	#endif
970
971	IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
972
973	u(1:mxcalc)=u(1:mxcalc) + timestep*( &
974	& du_phys(1:mxcalc) &
975	& +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) &
976	& +d_u_nudge(1:mxcalc) )
977	v(1:mxcalc)=v(1:mxcalc) + timestep*( &
978	& dv_phys(1:mxcalc) &
979	& +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) &
980	& +d_v_nudge(1:mxcalc) )
981	q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( &
982	& dq(1:mxcalc,:) &
983	& +d_q_adv(1:mxcalc,:) &
984	& +d_q_nudge(1:mxcalc,:) )
985
986	if (prt_level.ge.3) then
987	print *, &
988	& 'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ', &
989	& temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
990	print* ,'dv_phys=',dv_phys
991	print* ,'d_v_age=',d_v_age
992	print* ,'d_v_adv=',d_v_adv
993	print* ,'d_v_nudge=',d_v_nudge
994	print*, v
995	print*, vg
996	endif
997
998	temp(1:mxcalc)=temp(1:mxcalc)+timestep*( &
999	& dt_phys(1:mxcalc) &
1000	& +d_t_adv(1:mxcalc) &
1001	& +d_t_nudge(1:mxcalc) &
1002	& +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
1003
1004
1005	!=======================================================================
1006	!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
1007	!=======================================================================
1008
1009	teta=temp(pzero/play)*rkappa
1010
1011	!---------------------------------------------------------------------
1012	! Nudge soil temperature if requested
1013	!---------------------------------------------------------------------
1014
1015	IF (nudge_tsoil .AND. .NOT. lastcall) THEN
1016	ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) &
1017	& -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
1018	ENDIF
1019
1020	!---------------------------------------------------------------------
1021	! Add large-scale tendencies (advection, etc) :
1022	!---------------------------------------------------------------------
1023
1024	!cc nrlmd
1025	!cc tmpvar=teta
1026	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
1027	!cc
1028	!cc teta(1:mxcalc)=tmpvar(1:mxcalc)
1029	!cc tmpvar(:)=q(:,1)
1030	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
1031	!cc q(1:mxcalc,1)=tmpvar(1:mxcalc)
1032	!cc tmpvar(:)=q(:,2)
1033	!cc call advect_vert(llm,omega,timestep,tmpvar,plev)
1034	!cc q(1:mxcalc,2)=tmpvar(1:mxcalc)
1035
1036	END IF ! end if tendency of tendency should be added
1037
1038	!---------------------------------------------------------------------
1039	! Air temperature :
1040	!---------------------------------------------------------------------
1041	if (lastcall) then
1042	print*,'Pas de temps final ',it
1043	call ju2ymds(daytime, an, mois, jour, heure)
1044	print*,'a la date : a m j h',an, mois, jour ,heure/3600.
1045	endif
1046
1047	! incremente day time
1048	daytime = daytime+1./day_step
1049	day = int(daytime+0.1/day_step)
1050	! time = max(daytime-day,0.0)
1051	!Al1&jyg: correction de bug
1052	!cc time = real(mod(it,day_step))/day_step
1053	time = time_ini/24.+real(mod(it,day_step))/day_step
1054	it=it+1
1055
1056	enddo
1057
1058	if (ecrit_slab_oc.ne.-1) close(97)
1059
1060	!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
1061	! ---------------------------------------------------------------------------
1062	call dyn1dredem("restart1dyn.nc", &
1063	& plev,play,phi,phis,presnivs, &
1064	& u,v,temp,q,omega2)
1065
1066	CALL abort_gcm ('lmdz1d ','The End ',0)
1067
1068	END SUBROUTINE scm

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

Download in other formats:

Original Format