Context Navigation

scm.f90 @ 5446

Last change on this file since 5446 was 5433, checked in by evignon, 4 days ago

modification et correction de la prise en compte de l'advection verticale dans le 1D.
Il y a 3 aspects:

on tient compte de l'advection verticale des 3 phases de l'eau
on corrige l'application de la tendance qui etait calculee avant l'appel de la physique mais appliquee apres
on utilise desormais un schema amont par defaut (le schema centre, mis sous flag, conduisait a des oscillations). C'est un peu fort. Implementation a venir d'un schema van leer

ATTENTION: on perd la convergence sur les cas 1D avec vitesse verticale grande echelle (RICO, SANDU ...)

File size: 39.9 KB

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

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

Download in other formats:

Original Format