Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_lmd_driver.F @ 3067

Last change on this file since 3067 was 2758, checked in by aslmd, 3 years ago
moved the new driver as the main one in INTERFACES. saved the old one in the origin folder.
File size: 27.6 KB

Rev	Line
[1234]	1	!*******************************************************************************
	2	! PURPOSE: LMD_driver is the WRF mediation layer routine that provides the
	3	! interface to LMD physics packages in the WRF model layer. For those familiar
	4	! with the LMD GCM, the aim of this driver is to do part of the job of the
	5	! calfis.F routine.
	6	!*******************************************************************************
	7	! AUTHOR: A. Spiga - January 2007
	8	!*******************************************************************************
	9	! UPDATES: - included all final updates for the paper - March 2008
	10	! - general cleaning of code and comments - October 2008
	11	! - additions for idealized cases - January 2009
	12	! - additions for new soil model in physics - January 2010
	13	! - unified module_lmd_driver: old, new phys and LES - February 2011
	14	! - a new paradigm to prepare nested simulations - April 2014
[1580]	15	! - adapted to new interface philosophy & other planets - August 2016
[1234]	16	!*******************************************************************************
	17	MODULE module_lmd_driver
	18	CONTAINS
	19	SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, &
	20	IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, &
	21	i_start,i_end,j_start,j_end,kts,kte,num_tiles, &
	22	DX,DY, &
	23	MSFT,MSFU,MSFV, &
	24	GMT,JULYR,JULDAY, &
	25	P8W,DZ8W,T8W,Z,HT, &
	26	U,V,W,TH,T,P,EXNER,RHO, &
	27	PTOP, &
	28	RADT,CUDT, &
	29	TSK,PSFC, &
	30	RTHBLTEN,RUBLTEN,RVBLTEN, &
	31	num_3d_s,SCALAR, &
[2123]	32	num_3d_m,moist, &
[1234]	33	MARS_MODE, &
[1584]	34	planet_type, &
[1236]	35	M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
	36	M_H2OICE,M_TSOIL,M_Q2,M_TSURF, &
	37	M_FLUXRAD,M_WSTAR,M_ISOIL,M_DSOIL,&
	38	M_Z0, CST_Z0, M_GW, &
[1234]	39	NUM_SOIL_LAYERS, &
[1236]	40	CST_AL, CST_TI, &
	41	isfflx, diff_opt, km_opt, &
[1234]	42	HISTORY_INTERVAL, &
[2018]	43	HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,&
	44	CLOUDFRAC,TOTCLOUDFRAC, &
	45	GRAIN,GSNOW,REEVAP,SURFRAIN,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
	46	FLUXSURF_LW,FLXGRD,LSCEZ,H2OICE_REFF,LATENT_HF,SWDOWNZ,&
[1236]	47	TAU_DUST,RDUST,QSURFDUST,&
	48	MTOT,ICETOT,VMR_ICE,TAU_ICE,RICE,&
	49	HFMAX,ZMAX,&
	50	USTM,HFX,&
[1234]	51	SLPX,SLPY,RESTART)
	52	! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
	53
	54
	55
	56
	57
	58	!==================================================================
	59	! USES
	60	!==================================================================
	61	USE module_model_constants
	62	USE module_wrf_error
[1590]	63	!!!!!!!! interface modules
[1755]	64	USE variables_mod !! to set variables
[1577]	65	USE update_inputs_physiq_mod !! to set inputs for physiq
	66	USE update_outputs_physiq_mod !! to get outputs from physiq
[1579]	67	USE iniphysiq_mod !! to get iniphysiq subroutine
[1580]	68	USE callphysiq_mod !! to call the LMD physics
[1590]	69	!!!!!!!! interface modules
[1234]	70
	71	!==================================================================
	72	IMPLICIT NONE
	73	!==================================================================
	74
	75	!==================================================================
	76	! VARIABLES
	77	!==================================================================
	78
[1584]	79	CHARACTER(len=10),INTENT(IN) :: planet_type
[1234]	80	! WRF Dimensions
	81	INTEGER, INTENT(IN ) :: &
	82	ids,ide,jds,jde,kds,kde, &
	83	ims,ime,jms,jme,kms,kme, &
	84	kts,kte,num_tiles, &
	85	NUM_SOIL_LAYERS
	86	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
	87	i_start,i_end,j_start,j_end
	88	! Scalars
	89	INTEGER, INTENT(IN ) :: JULDAY, itimestep, julyr,id,max_dom,MARS_MODE
	90	INTEGER, INTENT(IN ) :: isfflx,diff_opt,km_opt
	91	REAL, INTENT(IN ) :: GMT,dt,dx,dy,RADT,CUDT
	92	REAL, INTENT(IN ) :: CST_AL, CST_TI
	93	REAL, INTENT(IN ) :: PTOP
[1236]	94	INTEGER, INTENT(IN ) :: HISTORY_INTERVAL
[1234]	95	! 2D arrays
[1236]	96	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
[1234]	97	MSFT,MSFU,MSFV, &
	98	XLAT,XLONG,HT, &
[1236]	99	M_ALBEDO,M_TI,M_EMISS, &
	100	SLPX,SLPY, &
	101	M_CO2ICE,M_H2OICE, &
	102	M_TSURF, M_Z0, &
	103	M_FLUXRAD,M_WSTAR, &
	104	PSFC,TSK
	105	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: &
	106	SWDOWNZ,&
	107	TAU_DUST,QSURFDUST,&
	108	MTOT,ICETOT,TAU_ICE,&
	109	HFMAX,ZMAX,&
[2018]	110	USTM,HFX,TOTCLOUDFRAC,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
	111	FLUXSURF_LW,FLXGRD,LATENT_HF,REEVAP,SURFRAIN
[1234]	112	! 3D arrays
	113	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
	114	dz8w,p8w,p,exner,t,t8w,rho,u,v,w,z,th
[1236]	115	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
	116	RTHBLTEN,RUBLTEN,RVBLTEN, &
[2018]	117	HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,RDUST,VMR_ICE,RICE,&
	118	CLOUDFRAC,GRAIN,GSNOW,LSCEZ,H2OICE_REFF
[1234]	119	REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
[1236]	120	M_Q2
[1234]	121	REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT ) :: &
[1236]	122	M_TSOIL,M_ISOIL, M_DSOIL
[1234]	123	REAL, INTENT(IN ) :: CST_Z0
	124	REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN ) :: &
[1236]	125	M_GW
[1234]	126	! 4D arrays
[2123]	127	INTEGER, INTENT(IN ) :: num_3d_s,num_3d_m
	128	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: scalar
	129	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_m ), INTENT(INOUT ) :: moist
[1234]	130	! Logical
	131	LOGICAL, INTENT(IN ) :: restart
	132
	133	!-------------------------------------------
	134	! LOCAL VARIABLES
	135	!-------------------------------------------
	136	INTEGER :: i,j,k,its,ite,jts,jte,ij,kk
	137	INTEGER :: subs,iii
	138
	139
[1580]	140	! * for tracer Mode 20 *
[1234]	141	REAL :: tau_decay
	142	! * ----------------------- *
	143
	144	! *** for LMD physics
	145	! ------> inputs:
	146	INTEGER :: ngrid,nlayer,nq,nsoil
[1724]	147	REAL :: MY
[1577]	148	REAL :: phisfi_val
[1578]	149	LOGICAL :: firstcall,lastcall
[1234]	150	! ----------
	151
	152	! <------ outputs:
	153	! physical tendencies
[1755]	154	REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdtheta
[1234]	155	! ... intermediate arrays
	156	REAL, DIMENSION(:), ALLOCATABLE :: &
	157	dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, &
	158	u_prof,v_prof,z_prof
	159	REAL, DIMENSION(:,:), ALLOCATABLE :: q_prof
	160
	161	! Additional control variables
	162	INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
[1724]	163	REAL :: elaps
[1234]	164	INTEGER :: test
	165	REAL :: wappel_phys
	166	LOGICAL, SAVE :: flag_first_restart
[1777]	167	INTEGER, SAVE :: previous_id = 0
[1234]	168	!**************************************************
	169	! IMPORTANT: pour travailler avec grid nesting
	170	! --- deux comportements distincts du save
	171	! ... ex1: ferme planeto avec PGF+MPI: standard
	172	! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE
	173	!**************************************************
	174	#ifdef SPECIAL_NEST_SAVE
	175	! une dimension supplementaire liee au nest
	176	REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
	177	dp_save
	178	REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
[1724]	179	du_save, dv_save, dt_save,dtheta_save
[1234]	180	REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
	181	dq_save
	182	#else
	183	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
	184	dp_save
	185	REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
[1724]	186	du_save, dv_save, dt_save,dtheta_save
[1234]	187	REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
	188	dq_save
	189	#endif
	190
	191	!!!IDEALIZED IDEALIZED
	192	REAL :: lat_input, lon_input, ls_input, lct_input
[1577]	193	LOGICAL :: isles
[1234]	194	!!!IDEALIZED IDEALIZED
	195
[1724]	196	REAL :: tk1,tk2
[1777]	197
[1234]	198	!==================================================================
	199	! CODE
	200	!==================================================================
	201
[1777]	202	print , '* ',planet_type,'** ENTER WRF-LMD DRIVER'
	203
	204	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1577]	205	!! determine here if this is turbulence-resolving or not
[1234]	206	IF (JULYR .ne. 9999) THEN
[1577]	207	isles = .false. ! "True" LES is not available in this version
[1234]	208	PRINT , ' REAL-CASE SIMULATION *'
	209	ELSE
	210	PRINT , ' IDEALIZED SIMULATION *'
	211	IF ((diff_opt .eq. 2) .and. (km_opt .eq. 2)) THEN
	212	PRINT , ' type: LES *'
	213	PRINT , ' diff_opt = 2 * km_opt = 2'
	214	PRINT , '** forcing is isfflx = ',isfflx
[1577]	215	isles = .true.
[1234]	216	!! SPECIAL LES
	217	ELSE
	218	PRINT , ' type: not LES *'
	219	PRINT , '** diff_opt = ',diff_opt
	220	PRINT , '** km_opt = ',km_opt
[1577]	221	isles = .false.
[1234]	222	!! IDEALIZED, no LES
	223	!! cependant, ne veut-on pas pouvoir
	224	!! prescrire un flux en idealise ??
	225	ENDIF
	226	ENDIF
[1777]	227	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1234]	228
	229	!-------------------------!
	230	! TWEAK on WRF DIMENSIONS !
	231	!-------------------------!
	232	its = i_start(1) ! define here one big tile
	233	ite = i_end(num_tiles)
	234	jts = j_start(1)
	235	jte = j_end(num_tiles)
	236	!!
[1577]	237	IF (isles .eqv. .false.) THEN
[1234]	238	relax=0
	239	sponge_top=0 ! another value than 0 triggers instabilities
	240	IF (id .gt. 1) relax=2 ! fix to avoid noise in nesting simulations ; 1 >> too much noise ...
	241	ENDIF
	242	ips=its
	243	ipe=ite
	244	jps=jts
	245	jpe=jte
[1577]	246	IF (isles .eqv. .false.) THEN
[1234]	247	IF (ips .eq. ids) ips=its+relax !! IF tests necesary for parallel runs
	248	IF (ipe .eq. ide-1) ipe=ite-relax
	249	IF (jps .eq. jds) jps=jts+relax
	250	IF (jpe .eq. jde-1) jpe=jte-relax
	251	ENDIF
	252	kps=kts !! start at surface
[1577]	253	IF (isles .eqv. .false.) THEN
[1234]	254	kpe=kte-sponge_top
	255	ELSE
	256	PRINT , ' IDEALIZED SIMULATION: LES * kpe=kte'
	257	kpe=kte !-sponge_top
	258	ENDIF
	259
	260	!----------------------------!
	261	! DIMENSIONS FOR THE PHYSICS !
	262	!----------------------------!
	263	!day_ini = JULDAY - 1 !! GCM convention !! pday at firstcall is day_ini
	264	wappel_phys = RADT
[1755]	265	zdt_split = dt*wappel_phys ! physical timestep (s)
[2743]	266	ptopwrf = ptop ! top pressure in WRF coordinates (in variables_mod)
[1234]	267	ngrid=(ipe-ips+1)*(jpe-jps+1) ! size of the horizontal grid
	268	nlayer = kpe-kps+1 ! number of vertical layers: nlayermx
	269	nsoil = NUM_SOIL_LAYERS ! number of soil layers: nsoilmx
	270	if (num_3d_s > 1) then ! number of advected fields
	271	nq = num_3d_s-1
	272	else
	273	nq = 1
	274	endif
	275	! **** needed but hardcoded
	276	lastcall = .false.
	277	! **** needed but hardcoded
	278
	279	!PRINT *, ips, ipe, jps, jpe
	280	!PRINT *, ngrid
	281
	282	elaps = (float(itimestep)-1.)*dt ! elapsed seconds of simulation
	283	!----------------------------------------------!
	284	! what is done at the first step of simulation !
	285	!----------------------------------------------!
	286	IF (elaps .eq. 0.) THEN
	287	flag_first_restart = .false.
	288	ELSE
	289	flag_first_restart=flag_first_restart.OR.(.NOT. ALLOCATED(dp_save))
	290	ENDIF
	291
[1777]	292	is_first_step: IF ((elaps .eq. 0.).or.flag_first_restart) THEN
[1234]	293	firstcall=.true. !! for continuity with GCM, physics are always called at the first WRF timestep
	294	!firstcall=.false. !! just in case you'd want to get rid of the physics
	295	test=0
	296	#ifdef SPECIAL_NEST_SAVE
	297	PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
	298	IF( .NOT. ALLOCATED( dp_save ) ) THEN
	299	PRINT , '* check ** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
	300	!! here are the arrays that would be useful to save physics tendencies
	301	ALLOCATE(dp_save(ngrid,max_dom))
	302	ALLOCATE(du_save(ngrid,nlayer,max_dom))
	303	ALLOCATE(dv_save(ngrid,nlayer,max_dom))
	304	ALLOCATE(dt_save(ngrid,nlayer,max_dom))
[1724]	305	ALLOCATE(dtheta_save(ngrid,nlayer,max_dom))
[1234]	306	ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
	307	dp_save(:,:)=0. !! initialize these arrays ...
	308	du_save(:,:,:)=0.
	309	dv_save(:,:,:)=0.
	310	dt_save(:,:,:)=0.
[1724]	311	dtheta_save(:,:,:)=0.
[1234]	312	dq_save(:,:,:,:)=0.
	313	ENDIF
	314	IF (id .lt. max_dom) THEN
	315	flag_first_restart=.true.
	316	ELSE
	317	flag_first_restart=.false.
	318	ENDIF
	319	#else
	320	IF( .NOT. ALLOCATED( dp_save ) ) THEN
	321	ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
	322	ALLOCATE(du_save(ngrid,nlayer))
	323	ALLOCATE(dv_save(ngrid,nlayer))
	324	ALLOCATE(dt_save(ngrid,nlayer))
[1724]	325	ALLOCATE(dtheta_save(ngrid,nlayer))
[1234]	326	ALLOCATE(dq_save(ngrid,nlayer,nq))
	327	ENDIF
	328	dp_save(:)=0. !! initialize these arrays ...
	329	du_save(:,:)=0.
	330	dv_save(:,:)=0.
	331	dt_save(:,:)=0.
[1724]	332	dtheta_save(:,:)=0.
[1234]	333	dq_save(:,:,:)=0.
	334	flag_first_restart=.false.
	335	#endif
	336	ELSE
	337	!-------------------------------------------------!
	338	! what is done for the other steps of simulation !
	339	!-------------------------------------------------!
	340	IF (wappel_phys .gt. 0.) THEN
	341	firstcall = .false.
	342	test = MODULO(itimestep-1,int(wappel_phys)) ! WRF time is integrated time (itimestep=1 at the end of first step)
	343	! -- same strategy as diagfi in the LMD GCM
	344	! -- arguments of modulo must be of the same kind (here: integers)
	345	ELSE
	346	firstcall = .false.
	347	test = 9999
	348	ENDIF
[1777]	349	ENDIF is_first_step
[1234]	350
[1777]	351	!------------------------------------!
	352	! print info about domain initially !
	353	! ... and whenever domain is changed !
	354	!------------------------------------!
	355	print ,'* ',planet_type,' ** DOMAIN',id
[1778]	356	IF (previous_id .ne. id) THEN
[1777]	357	print , '* ',planet_type,' ** ... INITIALIZE DOMAIN',id
[1778]	358	print , '* ',planet_type,' ** ... PREVIOUS DOMAIN was',previous_id
[1777]	359	print *, 'ITIMESTEP: ',itimestep
	360	print *, 'TILES: ', i_start,i_end, j_start, j_end ! numbers for simple runs, arrays for parallel runs
	361	print *, 'DOMAIN: ', ids, ide, jds, jde
	362	print *, 'MEMORY: ', ims, ime, jms, jme
	363	print *, 'COMPUT: ', ips, ipe, jps, jpe
	364	print *, 'SIZE OF PHYSICS GRID for this process: ',ngrid
	365	print *, 'ADVECTED TRACERS: ', num_3d_s-1
	366	print *, 'PHYSICS IS CALLED EACH...',wappel_phys
	367	print *, '-- means: PHYSICAL TIMESTEP=',zdt_split
	368	ENDIF
	369
	370
[1577]	371	!!!***********!!
	372	!!! IDEALIZED !!
	373	!!!***********!!
	374	IF (JULYR .eq. 9999) THEN
[1580]	375	PRINT ,'* ',planet_type,'** IDEALIZED SIMULATION'
	376	PRINT ,'* ',planet_type,'** BEWARE: input_coord must be here'
[1234]	377	open(unit=14,file='input_coord',form='formatted',status='old')
	378	rewind(14)
	379	read(14,*) lon_input
	380	read(14,*) lat_input
	381	read(14,*) ls_input
	382	read(14,*) lct_input
	383	close(14)
	384	ENDIF
	385
	386	!----------!
	387	! ALLOCATE !
	388	!----------!
[1755]	389	CALL allocate_interface(ngrid,nlayer,nq)
[1759]	390	ALLOCATE(pdtheta(ngrid,nlayer))
[1234]	391	!!!
	392	!!! BIG LOOP : 1. no call for physics, used saved values
	393	!!!
	394	IF (test.NE.0) THEN
[1580]	395	print ,'* ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test
[1234]	396	#ifdef SPECIAL_NEST_SAVE
[1755]	397	zdpsrf_omp(:)=dp_save(:,id)
	398	zdufi_omp(:,:)=du_save(:,:,id)
	399	zdvfi_omp(:,:)=dv_save(:,:,id)
	400	zdtfi_omp(:,:)=dt_save(:,:,id)
[1724]	401	pdtheta(:,:)=dtheta_save(:,:,id)
[1755]	402	zdqfi_omp(:,:,:)=dq_save(:,:,:,id)
[1234]	403	#else
[1724]	404	print*,'else'
[1755]	405	zdpsrf_omp(:)=dp_save(:)
	406	zdufi_omp(:,:)=du_save(:,:)
	407	zdvfi_omp(:,:)=dv_save(:,:)
	408	zdtfi_omp(:,:)=dt_save(:,:)
[1724]	409	pdtheta(:,:)=dtheta_save(:,:)
[1755]	410	zdqfi_omp(:,:,:)=dq_save(:,:,:)
[1234]	411	#endif
	412	!!!
	413	!!! BIG LOOP : 2. calculate physical tendencies
	414	!!!
	415	ELSE
	416	!----------!
	417	! ALLOCATE !
	418	!----------!
	419	! interm
	420	ALLOCATE(dz8w_prof(nlayer))
	421	ALLOCATE(p8w_prof(nlayer))
	422	ALLOCATE(p_prof(nlayer))
	423	ALLOCATE(t_prof(nlayer))
	424	ALLOCATE(t8w_prof(nlayer))
	425	ALLOCATE(u_prof(nlayer))
	426	ALLOCATE(v_prof(nlayer))
	427	ALLOCATE(z_prof(nlayer))
	428	ALLOCATE(q_prof(nlayer,nq))
	429
	430
	431	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	432	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	433	!! PREPARE PHYSICS INPUTS !!
	434	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	435	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	436
[1577]	437	!! INITIALIZE AND ALLOCATE EVERYTHING !! here, only firstcall
[1234]	438	allocation_firstcall: IF (firstcall .EQV. .true.) THEN
[1577]	439	!! tracers' name
[1580]	440	PRINT ,'* ',planet_type,'** TRACERS NAMES'
[1740]	441	CALL update_inputs_physiq_tracers(nq,MARS_MODE)
[1234]	442	!! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM)
	443	!! parameters are defined in the module_model_constants.F WRF routine
[1580]	444	PRINT ,'* ',planet_type,'** INITIALIZE ARRAYS FOR PHYSICS'
[1578]	445	!! need to get initial time first
	446	CALL update_inputs_physiq_time(&
	447	JULYR,JULDAY,GMT,&
	448	elaps,&
	449	lct_input,lon_input,ls_input,&
[1755]	450	MY)
[1759]	451	!! Set up initial time
	452	phour_ini = JH_cur_split
[1234]	453	!! Fill planetary parameters in modules
	454	!! Values defined in the module_model_constants.F WRF routine
[1577]	455	CALL update_inputs_physiq_constants
[1579]	456	!! Initialize physics
[1755]	457	CALL iniphysiq(ngrid,nlayer,nq,wappel_phys,&
	458	wdaysec,floor(JD_cur), &
[1579]	459	1./reradius,g,r_d,cp,1)
[2070]	460
	461
[1577]	462	ENDIF allocation_firstcall
[1234]	463
	464	!!*****************************!!
	465	!! PREPARE "FOOD" FOR PHYSIQ.F !!
	466	!!*****************************!!
	467
	468	DO j = jps,jpe
	469	DO i = ips,ipe
	470
	471	!!*******************************************!!
	472	!! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !!
	473	!!*******************************************!!
	474
	475	!-----------------------------------!
	476	! 1D subscript for physics "cursor" !
	477	!-----------------------------------!
	478	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
	479
	480	!--------------------------------------!
	481	! 1D columns : inputs for the physics !
	482	! ... vert. coord., meteor. fields ... !
	483	!--------------------------------------!
	484	dz8w_prof(:) = dz8w(i,kps:kpe,j) ! dz between full levels (m)
[1755]	485	p8w_prof(:) = p8w(i,kps:kpe,j) ! pressure full level (Pa) >> zplev_omp
	486	p_prof(:) = p(i,kps:kpe,j) ! pressure half level (Pa) >> zplay_omp
[1234]	487	t_prof(:) = t(i,kps:kpe,j) ! temperature half level (K) >> pt
	488	t8w_prof(:) = t8w(i,kps:kpe,j) ! temperature full level (K)
[1755]	489	u_prof(:) = u(i,kps:kpe,j) ! zonal wind (A-grid: unstaggered) half level (m/s) >> zufi_omp
[1234]	490	v_prof(:) = v(i,kps:kpe,j) ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
[1755]	491	z_prof(:) = z(i,kps:kpe,j) ! geopotential height half level (m) >> zphi_omp/g
[1234]	492
	493	!--------------------------------!
	494	! specific treatment for tracers !
	495	!--------------------------------!
	496	IF (MARS_MODE .EQ. 0) THEN
	497	q_prof(:,1)=0.95
[2208]	498	ELSE IF (MARS_MODE .EQ. 42) THEN
[2123]	499	q_prof(:,1) = moist(i,kps:kpe,j,P_QV) !! the names were set above !! one dummy tracer in WRF
	500	q_prof(:,2) = SCALAR(i,kps:kpe,j,3) !! the names were set above !! one dummy tracer in WRF
[2208]	501	ELSE
	502	q_prof(:,1:nq) = SCALAR(i,kps:kpe,j,2:nq+1) !! the names were set above !! one dummy tracer in WRF
[1234]	503	ENDIF
	504
	505	IF (MARS_MODE .EQ. 20) THEN
	506	IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
	507	q_prof(:,:) = 0.95
	508	ENDIF
	509	ENDIF
	510
[1741]	511	IF (MARS_MODE .EQ. 32) THEN
	512	IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
	513	q_prof(:,7) = 0.95
	514	!! traceurs(7) = 'co2'
	515	ENDIF
	516	ENDIF
	517
[1234]	518	IF (firstcall .EQV. .true.) THEN
	519	!-----------------!
	520	! Optional output !
	521	!-----------------!
	522	IF ( (i == ips) .AND. (j == jps) ) THEN
	523	PRINT *,'z_prof ',z_prof
	524	PRINT *,'dz8w_prof',dz8w_prof
	525	PRINT *,'p8w_prof ',p8w_prof
	526	PRINT *,'p_prof ',p_prof
	527	PRINT *,'t_prof ',t_prof
	528	PRINT *,'t8w_prof ',t8w_prof
	529	PRINT *,'u_prof ',u_prof
	530	PRINT *,'v_prof ',v_prof
	531	PRINT *,'q_prof ',q_prof
	532	ENDIF
	533	ENDIF
	534
	535	!---------------------------------------------!
	536	! in LMD physics, geopotential must be !
	537	! expressed with respect to the local surface !
	538	!---------------------------------------------!
[1755]	539	zphi_omp(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) )
[1234]	540
	541	!--------------------------------!
	542	! Dynamic fields for LMD physics !
	543	!--------------------------------!
[1755]	544	zplev_omp(subs,1:nlayer) = p8w_prof(1:nlayer) !! NB: last level: no data
	545	zplay_omp(subs,:) = p_prof(:)
	546	ztfi_omp(subs,:) = t_prof(:)
	547	zufi_omp(subs,:) = u_prof(:)
	548	zvfi_omp(subs,:) = v_prof(:)
	549	flxwfi_omp(subs,:) = 0 !! NB: not used in the physics, only diagnostic...
[1234]	550	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	551	!! for IDEALIZED CASES ONLY
[1755]	552	IF (JULYR .eq. 9999) zplev_omp(subs,nlayer+1)=0. !! zplev_omp(subs,nlayer+1)=ptop >> NO !
[1234]	553	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	554
	555	! NOTE:
[1755]	556	! IF ( zplev_omp(subs,nlayer+1) .le. 0 ) zplev_omp(subs,nlayer+1)=ptop
[1234]	557	! cree des diagnostics delirants et aleatoires dans le transfert radiatif
	558
	559	!---------!
	560	! Tracers !
	561	!---------!
[1755]	562	zqfi_omp(subs,:,:) = q_prof(:,:) !! traceurs generiques, seuls noms sont specifiques
[1234]	563
	564	ENDDO
	565	ENDDO
	566
	567	!---------------------------------------------------------!
	568	! Ground geopotential !
	569	! (=g*HT(i,j), only used in the microphysics: newcondens) !
	570	!---------------------------------------------------------!
	571	phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.) !! NB: z_prof are half levels, so z_prof(1) is not the surface
	572
	573	!!*****************!!
	574	!! CLEAN THE PLACE !!
	575	!!*****************!!
	576	DEALLOCATE(dz8w_prof)
	577	DEALLOCATE(z_prof)
	578	DEALLOCATE(p8w_prof)
	579	DEALLOCATE(p_prof)
	580	DEALLOCATE(t_prof)
	581	DEALLOCATE(t8w_prof)
	582	DEALLOCATE(u_prof)
	583	DEALLOCATE(v_prof)
	584	DEALLOCATE(q_prof)
	585
[1577]	586
	587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	588	!!! ONE DOMAIN CASE
	589	!!! --> we simply need to initialize static and physics inputs
	590	!!! SEVERAL DOMAINS
	591	!!! --> we update static and physics inputs each time
	592	!!! to account for domain change
	593	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	594	pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN
[1580]	595	PRINT ,'* ',planet_type,'** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.'
[1577]	596	!!*******************************************!!
	597	!!*******************************************!!
	598	CALL update_inputs_physiq_geom( &
	599	ims,ime,jms,jme,&
	600	ips,ipe,jps,jpe,&
[1595]	601	JULYR,ngrid,nlayer,&
[1577]	602	DX,DY,MSFT,&
	603	lat_input, lon_input,&
	604	XLAT,XLONG)
	605	!!!
	606	CALL update_inputs_physiq_slope( &
	607	ims,ime,jms,jme,&
	608	ips,ipe,jps,jpe,&
	609	JULYR,&
	610	SLPX,SLPY)
	611	!!!
	612	CALL update_inputs_physiq_soil( &
	613	ims,ime,jms,jme,&
	614	ips,ipe,jps,jpe,&
	615	JULYR,nsoil,&
	616	M_TI,CST_TI,&
	617	M_ISOIL,M_DSOIL,&
	618	M_TSOIL,M_TSURF)
	619	!!!
	620	CALL update_inputs_physiq_surf( &
	621	ims,ime,jms,jme,&
	622	ips,ipe,jps,jpe,&
	623	JULYR,MARS_MODE,&
	624	M_ALBEDO,CST_AL,&
	625	M_TSURF,M_EMISS,M_CO2ICE,&
[1742]	626	M_GW,M_Z0,CST_Z0,&
[1577]	627	M_H2OICE,&
	628	phisfi_val)
	629	!!!
	630	CALL update_inputs_physiq_turb( &
	631	ims,ime,jms,jme,kms,kme,&
	632	ips,ipe,jps,jpe,&
	633	RESTART,isles,&
	634	M_Q2,M_WSTAR)
	635	!!!
	636	CALL update_inputs_physiq_rad( &
	637	ims,ime,jms,jme,&
	638	ips,ipe,jps,jpe,&
	639	RESTART,&
	640	M_FLUXRAD)
	641	!!!
	642	ENDIF pass_interface
	643	!!*******************************************!!
	644	!!*******************************************!!
	645
[1234]	646	!!!!!!!!!!!!!!!!!!!!!!
	647	!!!!!!!!!!!!!!!!!!!!!!
	648	!! CALL LMD PHYSICS !!
	649	!!!!!!!!!!!!!!!!!!!!!!
	650	!!!!!!!!!!!!!!!!!!!!!!
	651
	652	call_physics : IF (wappel_phys .ne. 0.) THEN
[1577]	653	!!! initialize tendencies (security)
[1755]	654	zdpsrf_omp(:)=0.
	655	zdufi_omp(:,:)=0.
	656	zdvfi_omp(:,:)=0.
	657	zdtfi_omp(:,:)=0.
[1724]	658	pdtheta(:,:)=0.
[1755]	659	zdqfi_omp(:,:,:)=0.
[1580]	660	print , '* ',planet_type,'** CALL TO LMD PHYSICS'
[1577]	661	!!!
	662	CALL update_inputs_physiq_time(&
	663	JULYR,JULDAY,GMT,&
	664	elaps,&
	665	lct_input,lon_input,ls_input,&
[1755]	666	MY)
[1577]	667	!!!
[1739]	668	CALL call_physiq(planet_type,ngrid,nlayer,nq, &
[1755]	669	firstcall,lastcall)
[1580]	670	!!!
[1411]	671
	672	!! specific scenario. necessary to add the right amount of dust.
	673	#ifdef DUSTSTORM
	674	IF (firstcall .EQV. .true.) THEN
[1755]	675	zdqfi_omp(:,:,:) = zdqfi_omp(:,:,:) / dt
[1411]	676	ENDIF
	677	#endif
	678
[1724]	679	IF (planet_type .eq. "venus" ) THEN
	680	DO j=jps,jpe
	681	DO i=ips,ipe
	682	do k=kps,kpe
	683	subs=(j-jps)*(ipe-ips+1)+(i-ips+1)
[1755]	684	tk1=(ztfi_omp(subs,k)*nu + nuTT00*nulog((p1000mb/zplay_omp(subs,k))rcp))(1/nu)
	685	tk2=((ztfi_omp(subs,k) + zdtfi_omp(subs,k))*nu + nuTT00*nulog((p1000mb/zplay_omp(subs,k))rcp))(1/nu)
[1724]	686	pdtheta(subs,k)=tk2-tk1
	687	enddo
	688	ENDDO
	689	ENDDO
	690	ENDIF
	691
[1580]	692	print , '* ',planet_type,'** CALL TO LMD PHYSICS DONE'
[1234]	693
	694	!---------------------------------------------------------------------------------!
	695	! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
	696	!---------------------------------------------------------------------------------!
	697	#ifdef SPECIAL_NEST_SAVE
[1755]	698	dp_save(:,id)=zdpsrf_omp(:)
	699	du_save(:,:,id)=zdufi_omp(:,:)
	700	dv_save(:,:,id)=zdvfi_omp(:,:)
	701	dt_save(:,:,id)=zdtfi_omp(:,:)
[1724]	702	dtheta_save(:,:,id)=pdtheta(:,:)
[1755]	703	dq_save(:,:,:,id)=zdqfi_omp(:,:,:)
[1234]	704	#else
[1755]	705	dp_save(:)=zdpsrf_omp(:)
	706	du_save(:,:)=zdufi_omp(:,:)
	707	dv_save(:,:)=zdvfi_omp(:,:)
	708	dt_save(:,:)=zdtfi_omp(:,:)
[1724]	709	dtheta_save(:,:)=pdtheta(:,:)
[1755]	710	dq_save(:,:,:)=zdqfi_omp(:,:,:)
[1234]	711	#endif
[1236]	712
	713	!! OUTPUT OUTPUT OUTPUT
[1577]	714	!-------------------------------------------------------!
	715	! Save key variables for restart and output and nesting !
	716	!-------------------------------------------------------!
	717	!!!
	718	CALL update_outputs_physiq_surf( &
	719	ims,ime,jms,jme,&
	720	ips,ipe,jps,jpe,&
	721	MARS_MODE,&
	722	M_TSURF,M_CO2ICE,&
	723	M_H2OICE)
	724	!!!
	725	CALL update_outputs_physiq_soil( &
	726	ims,ime,jms,jme,&
	727	ips,ipe,jps,jpe,&
	728	nsoil,&
	729	M_TSOIL)
	730	!!!
	731	CALL update_outputs_physiq_rad( &
	732	ims,ime,jms,jme,&
	733	ips,ipe,jps,jpe,&
	734	M_FLUXRAD)
	735	!!!
	736	CALL update_outputs_physiq_turb( &
	737	ims,ime,jms,jme,kms,kme,&
[1724]	738	ips,ipe,jps,jpe,kps,kpe,&
[1577]	739	M_Q2,M_WSTAR,&
	740	HFMAX,ZMAX,USTM,HFX)
[1590]	741	!!!
	742	CALL update_outputs_physiq_diag( &
	743	ims,ime,jms,jme,kms,kme,&
	744	ips,ipe,jps,jpe,kps,kpe,&
	745	SWDOWNZ,TAU_DUST,QSURFDUST,&
	746	MTOT,ICETOT,TAU_ICE,&
[2018]	747	HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,&
	748	RDUST,VMR_ICE,RICE,&
	749	CLOUDFRAC,TOTCLOUDFRAC,&
	750	GRAIN,GSNOW,REEVAP,SURFRAIN,&
	751	ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
	752	FLUXSURF_LW,FLXGRD,LSCEZ,H2OICE_REFF,LATENT_HF)
[1590]	753	!!!
[1777]	754	print , '* ',planet_type,'** OUTPUT PHYSICS DONE'
[1724]	755
[1234]	756	ENDIF call_physics
	757
	758	ENDIF ! end of BIG LOOP 2.
	759	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	760	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	761	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	762	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	763
	764	!!***************************!!
	765	!! DEDUCE TENDENCIES FOR WRF !!
	766	!!***************************!!
	767	RTHBLTEN(ims:ime,kms:kme,jms:jme)=0.
	768	RUBLTEN(ims:ime,kms:kme,jms:jme)=0.
	769	RVBLTEN(ims:ime,kms:kme,jms:jme)=0.
	770	PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF
	771	!------------------------------------------------------------------!
	772	! WRF adds the physical and dynamical tendencies !
	773	! thus the tendencies must be passed as 'per second' tendencies !
	774	! --when physics is not called, the applied physical tendency !
	775	! --is the one calculated during the last call to physics !
	776	!------------------------------------------------------------------!
[1724]	777	!print*,'pdt',pdt(1,1),pdt(1,nlayer)
	778	!print*,'exner',exner(1,:,1)
[1234]	779	DO j = jps,jpe
	780	DO i = ips,ipe
	781
	782	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
	783
	784	! zonal wind
[1755]	785	RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe)
[1234]	786	! meridional wind
[1755]	787	RVBLTEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe)
[1234]	788	! potential temperature
	789	! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
[1724]	790	IF (planet_type .eq. "venus" ) THEN
	791	RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe)
	792	ELSE
[1759]	793	RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
[1724]	794	ENDIF
[1234]	795	! update surface pressure (cf CO2 cycle in physics)
	796	! here dt is needed
[1755]	797	PSFC(i,j)=PSFC(i,j)+zdpsrf_omp(subs)*dt
[1234]	798	! tracers
	799	SCALAR(i,kps:kpe,j,1)=0.
	800	SELECT CASE (MARS_MODE)
	801	CASE(0)
	802	SCALAR(i,kps:kpe,j,:)=0.
	803	CASE(20)
[1580]	804	!! tracer mode 20 : add a passive tracer with radioactive-like decay
[1234]	805	IF ( (i == ips) .AND. (j == jps) ) print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER.'
	806	tau_decay=60.*10. !! why not make it a namelist argument?
	807	SCALAR(i,kps:kpe,j,2) = SCALAR(i,kps:kpe,j,2)*exp(-dt/tau_decay)
	808	SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer
[2208]	809	CASE(42)
[2123]	810	moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV)+zdqfi_omp(subs,kps:kpe,1)*dt !!! here dt is needed
	811	scalar(i,kps:kpe,j,3)=scalar(i,kps:kpe,j,3)+zdqfi_omp(subs,kps:kpe,2)*dt !!! here dt is needed
[2208]	812	CASE DEFAULT
	813	SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+zdqfi_omp(subs,kps:kpe,1:nq)*dt !!! here dt is needed
[1234]	814	END SELECT
[2208]	815
[1234]	816	ENDDO
	817	ENDDO
[1755]	818	CALL deallocate_interface
[1724]	819	DEALLOCATE(pdtheta)
[1234]	820	!!*****!!
	821	!! END !!
	822	!!*****!!
[1580]	823	print ,'* ',planet_type,'** END LMD PHYSICS'
[1777]	824	previous_id = id
[1234]	825	!----------------------------------------------------------------!
	826	! use debug (see solve_em) if you wanna display some message ... !
	827	!----------------------------------------------------------------!
	828	END SUBROUTINE lmd_driver
	829
	830	END MODULE module_lmd_driver

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

Download in other formats:

Original Format