Context Navigation

lmdz1d.F90 @ 2053

Last change on this file since 2053 was 2040, checked in by Ehouarn Millour, 11 years ago

Added a "strato_custom" mode in disvert (for development purposes).
This is triggered by setting "vert_sampling=strato_custom" in a def file.

More work is needed to clean disvert (and merge it with disvert_noterre).
DC+EM

File size: 36.9 KB

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

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

Download in other formats:

Original Format