Context Navigation

scm.f90 @ 5623

Last change on this file since 5623 was 5609, checked in by aborella, 3 months ago

changed treatment of prognostic variables for prognostic clouds
adapted sedimentation and autoconversion for prognostic cirrus clouds
cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
simplified assumptions for ice crystals deposition / sublimation
first version of the coupling between prognostic cirrus clouds and deep convection
added persistent contrail cirrus clouds in radiative diagnostics

File size: 40.5 KB

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

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

Download in other formats:

Original Format