Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

update_inputs_physiq_mod.F

Last change on this file was 1755, checked in by aslmd, 7 years ago

MESOSCALE MARS and VENUS. moved the interface arrays in a module named variables_mod.F in dynphy_wrf to get a greater flexibility for the various planets (regarding e.g. REAL 4 or 8). incidentally this makes module_lmd_driver being more modular and the interfacing between dynamics and physics is further simplified by the fields shared in variables_mod. compilation was checked on Mars but not on Venus

File size: 23.7 KB

Line
1	MODULE update_inputs_physiq_mod
2
3	IMPLICIT NONE
4
5	CHARACTER(len=20),save,allocatable,dimension(:) :: traceurs ! tracer names
6
7	CONTAINS
8
9	!SUBROUTINE update_inputs_physiq_time
10	!SUBROUTINE update_inputs_physiq_tracers
11	!SUBROUTINE update_inputs_physiq_constants
12	!SUBROUTINE update_inputs_physiq_geom
13	!SUBROUTINE update_inputs_physiq_surf
14	!SUBROUTINE update_inputs_physiq_soil
15	!SUBROUTINE update_inputs_physiq_turb
16	!SUBROUTINE update_inputs_physiq_rad
17	!SUBROUTINE update_inputs_physiq_slope
18
19	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
20	real function ls2sol(ls)
21
22	!c Returns solar longitude, Ls (in deg.), from day number (in sol),
23	!c where sol=0=Ls=0 at the northern hemisphere spring equinox
24
25	implicit none
26
27	!c Arguments:
28	real ls
29
30	!c Local:
31	double precision xref,zx0,zteta,zz
32	!c xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
33	double precision year_day
34	double precision peri_day,timeperi,e_elips
35	double precision pi,degrad
36	parameter (year_day=668.6d0) ! number of sols in a amartian year
37	!c data peri_day /485.0/
38	parameter (peri_day=485.35d0) ! date (in sols) of perihelion
39	!c timeperi: 2pi( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
40	parameter (timeperi=1.90258341759902d0)
41	parameter (e_elips=0.0934d0) ! eccentricity of orbit
42	parameter (pi=3.14159265358979d0)
43	parameter (degrad=57.2957795130823d0)
44
45	if (abs(ls).lt.1.0e-5) then
46	if (ls.ge.0.0) then
47	ls2sol = 0.0
48	else
49	ls2sol = year_day
50	end if
51	return
52	end if
53
54	zteta = ls/degrad + timeperi
55	zx0 = 2.0datan(dtan(0.5zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
56	xref = zx0-e_elips*dsin(zx0)
57	zz = xref/(2.*pi)
58	ls2sol = zz*year_day + peri_day
59	if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day
60	if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day
61
62	return
63	end function ls2sol
64	!!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
65
66	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
67	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
68	SUBROUTINE update_inputs_physiq_time(&
69	JULYR,JULDAY,GMT,&
70	elaps,&
71	lct_input,lon_input,ls_input,&
72	MY)
73
74	USE variables_mod, only: JD_cur,JH_cur_split,phour_ini
75	!! JD_cur <> pday ! Julian day
76	!! JH_cur_split <> ptime ! Julian hour (fraction of day)
77
78	implicit none
79
80	INTEGER, INTENT(IN) :: JULDAY, JULYR
81	REAL, INTENT(IN) :: GMT,elaps,lon_input,ls_input,lct_input
82	REAL,INTENT(OUT) :: MY
83
84	IF (JULYR .ne. 9999) THEN
85	!
86	! specified
87	!
88	JH_cur_split = (GMT + elaps/3700.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
89	JH_cur_split = MODULO(JH_cur_split,24.) !! the two arguments of MODULO must be of the same type
90	JH_cur_split = JH_cur_split / 24.
91	JD_cur = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
92	JD_cur = MODULO(int(JD_cur),669)
93	MY = (JULYR-2000) + (88800.(JULDAY - 1)+3700.GMT+elaps)/59496000.
94	MY = INT(MY)
95	ELSE
96	!
97	! idealized
98	!
99	JH_cur_split = lct_input - lon_input / 15. + elaps/3700.
100	JH_cur_split = MODULO(JH_cur_split,24.)
101	JH_cur_split = JH_cur_split / 24.
102	JD_cur = floor(ls2sol(ls_input)) + INT((3700*(lct_input - lon_input / 15.) + elaps)/88800)
103	JD_cur = MODULO(int(JD_cur),669)
104	MY = 2024
105	!day_ini = floor(ls2sol(ls_input)) !! JD_cur at firstcall is day_ini
106	ENDIF
107	print ,'* Mars ** TIME IS', JD_cur, JH_cur_split*24.
108
109	END SUBROUTINE update_inputs_physiq_time
110
111	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
112	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
113	SUBROUTINE update_inputs_physiq_tracers(nq,MARS_MODE)
114
115	INTEGER, INTENT(IN) :: nq,MARS_MODE
116
117	ALLOCATE(traceurs(nq))
118
119	!!! name of tracers to mimic entries in tracer.def
120	!!! ----> IT IS IMPORTANT TO KEEP SAME ORDERING AS IN THE REGISTRY !!!!
121	!!! SAME NAMING AS IN THE LMD PHYSICS !!!!
122	SELECT CASE (MARS_MODE)
123	CASE(0,10)
124	traceurs(nq) = 'co2'
125	CASE(1) ! scalar:qh2o,qh2o_ice
126	traceurs(1) = 'h2o_vap'
127	traceurs(2) = 'h2o_ice'
128	CASE(2) ! scalar:qdust
129	traceurs(1) = 'dust01'
130	CASE(3) ! scalar:qdust,qdustn
131	traceurs(1) = 'dust_mass'
132	traceurs(2) = 'dust_number'
133	CASE(11) ! scalar:qh2o,qh2o_ice,qdust,qdustn
134	traceurs(1) = 'h2o_vap'
135	traceurs(2) = 'h2o_ice'
136	traceurs(3) = 'dust_mass'
137	traceurs(4) = 'dust_number'
138	CASE(12)
139	traceurs(1) = 'h2o_vap'
140	traceurs(2) = 'h2o_ice'
141	traceurs(3) = 'dust_mass'
142	traceurs(4) = 'dust_number'
143	traceurs(5) = 'ccn_mass'
144	traceurs(6) = 'ccn_number'
145	CASE(20)
146	traceurs(1) = 'qtrac1'
147	CASE(32) ! scalar:qh2o,qh2o_ice,qdust,qdustn,qccn,qccnn,qco2,qco2_ice,qccn_co2,qccnn_co2
148	traceurs(1) = 'h2o_vap'
149	traceurs(2) = 'h2o_ice'
150	traceurs(3) = 'dust_mass'
151	traceurs(4) = 'dust_number'
152	traceurs(5) = 'ccn_mass'
153	traceurs(6) = 'ccn_number'
154	traceurs(7) = 'co2'
155	traceurs(8) = 'co2_ice'
156	traceurs(9) = 'ccnco2_mass'
157	traceurs(10) = 'ccnco2_number'
158	CASE(42)
159	traceurs(1) = 'co2'
160	traceurs(2) = 'co'
161	traceurs(3) = 'o'
162	traceurs(4) = 'o1d'
163	traceurs(5) = 'o2'
164	traceurs(6) = 'o3'
165	traceurs(7) = 'h'
166	traceurs(8) = 'h2'
167	traceurs(9) = 'oh'
168	traceurs(10) = 'ho2'
169	traceurs(11) = 'h2o2'
170	traceurs(12) = 'ch4'
171	traceurs(13) = 'n2'
172	traceurs(14) = 'ar'
173	traceurs(15) = 'h2o_ice'
174	traceurs(16) = 'h2o_vap'
175	traceurs(17) = 'dust_mass'
176	traceurs(18) = 'dust_number'
177	END SELECT
178
179	!!---------------------!!
180	!! OUTPUT FOR CHECKING !!
181	!!---------------------!!
182	print*,"check: traceurs",traceurs
183
184	END SUBROUTINE update_inputs_physiq_tracers
185
186	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188	SUBROUTINE update_inputs_physiq_constants
189
190	USE module_model_constants
191	use comcstfi_h, only: omeg,mugaz
192	use planete_h, only: year_day,periheli,aphelie, &
193	peri_day,obliquit,emin_turb, &
194	lmixmin
195	use surfdat_h, only: emissiv,albedice,iceradius, &
196	emisice,dtemisice, &
197	z0_default
198	use comsoil_h, only: volcapa
199
200	!! comcstfi_h
201	omeg = womeg
202	mugaz = wmugaz
203	print*,"check: omeg,mugaz"
204	print*,omeg,mugaz
205	!! planet_h
206	year_day = wyear_day
207	periheli = wperiheli
208	aphelie = waphelie
209	peri_day = wperi_day
210	obliquit = wobliquit
211	emin_turb = wemin_turb
212	lmixmin = wlmixmin
213	print*,"check: year_day,periheli,aphelie,peri_day,obliquit"
214	print*,year_day,periheli,aphelie,peri_day,obliquit
215	print*,"check: emin_turb,lmixmin"
216	print*,emin_turb,lmixmin
217	!! surfdat_h
218	emissiv = wemissiv
219	emisice(1) = wemissiceN
220	emisice(2) = wemissiceS
221	albedice(1) = walbediceN
222	albedice(2) = walbediceS
223	iceradius(1) = wiceradiusN
224	iceradius(2) = wiceradiusS
225	dtemisice(1) = wdtemisiceN
226	dtemisice(2) = wdtemisiceS
227	z0_default = wz0
228	print*,"check: z0def,emissiv,emisice,albedice,iceradius,dtemisice"
229	print*,z0_default,emissiv,emisice,albedice,iceradius,dtemisice
230	!! comsoil_h
231	volcapa = wvolcapa
232	print*,"check: volcapa",volcapa
233
234	END SUBROUTINE update_inputs_physiq_constants
235
236	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
237	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238	SUBROUTINE update_inputs_physiq_geom( &
239	ims,ime,jms,jme,&
240	ips,ipe,jps,jpe,&
241	JULYR,ngrid,nlayer,&
242	DX,DY,MSFT,&
243	lat_input, lon_input,&
244	XLAT,XLONG)
245
246	! in WRF (share)
247	USE module_model_constants, only: DEGRAD
248	! in LMD (phymars)
249	use comgeomfi_h, only: ini_fillgeom
250	! in LMD (phy_common)
251	USE mod_grid_phy_lmdz, ONLY: init_grid_phy_lmdz
252	USE geometry_mod, ONLY: latitude,latitude_deg,&
253	longitude,longitude_deg,&
254	cell_area
255
256	INTEGER, INTENT(IN) :: ims,ime,jms,jme
257	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,ngrid,nlayer
258	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
259	MSFT,XLAT,XLONG
260	REAL, INTENT(IN) :: dx,dy
261	REAL, INTENT(IN) :: lat_input, lon_input
262	INTEGER :: i,j,subs
263	REAL,DIMENSION(ngrid) :: plon, plat, parea
264
265	DO j = jps,jpe
266	DO i = ips,ipe
267
268	!-----------------------------------!
269	! 1D subscript for physics "cursor" !
270	!-----------------------------------!
271	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
272
273	!----------------------------------------!
274	! Surface of each part of the grid (m^2) !
275	!----------------------------------------!
276	!parea(subs) = dx*dy !! 1. idealized cases - computational grid
277	parea(subs) = (dx/msft(i,j))*(dy/msft(i,j)) !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
278	!parea(subs)=dx*dy/msfu(i,j) !! 3. special for Mercator GCM-like simulations
279
280	!---------------------------------------------!
281	! Mass-point latitude and longitude (radians) !
282	!---------------------------------------------!
283	IF (JULYR .ne. 9999) THEN
284	plat(subs) = XLAT(i,j)*DEGRAD
285	plon(subs) = XLONG(i,j)*DEGRAD
286	ELSE
287	!!! IDEALIZED CASE
288	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION lat: ',lat_input
289	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION lon: ',lon_input
290	plat(subs) = lat_input*DEGRAD
291	plon(subs) = lon_input*DEGRAD
292	ENDIF
293
294	ENDDO
295	ENDDO
296
297	!! FILL GEOMETRICAL ARRAYS !!
298	call ini_fillgeom(ngrid,plat,plon,parea)
299
300	!!! ----------------------------------------------------------
301	!!! --- initializing geometry in phy_common
302	!!! --- (this is quite planet-independent)
303	!!! ----------------------------------------------------------
304	! initialize mod_grid_phy_lmdz
305	CALL init_grid_phy_lmdz(1,1,ipe-ips+1,jpe-jps+1,nlayer)
306	! fill in geometry_mod variables
307	! ... copy over local grid longitudes and latitudes
308	! ... partly what is done in init_geometry
309	IF(.not.ALLOCATED(longitude)) ALLOCATE(longitude(ngrid))
310	IF(.not.ALLOCATED(longitude_deg)) ALLOCATE(longitude_deg(ngrid))
311	IF(.not.ALLOCATED(latitude)) ALLOCATE(latitude(ngrid))
312	IF(.not.ALLOCATED(latitude_deg)) ALLOCATE(latitude_deg(ngrid))
313	IF(.not.ALLOCATED(cell_area)) ALLOCATE(cell_area(ngrid))
314	longitude(:) = plon(:)
315	latitude(:) = plat(:)
316	longitude_deg(:) = plon(:)/DEGRAD
317	latitude_deg(:) = plat(:)/DEGRAD
318	cell_area(:) = parea(:)
319	!!! ----------------------------------------------------------
320
321	END SUBROUTINE update_inputs_physiq_geom
322
323	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
325	SUBROUTINE update_inputs_physiq_surf( &
326	ims,ime,jms,jme,&
327	ips,ipe,jps,jpe,&
328	JULYR,MARS_MODE,&
329	M_ALBEDO,CST_AL,&
330	M_TSURF,M_EMISS,M_CO2ICE,&
331	M_GW,M_Z0,CST_Z0,&
332	M_H2OICE,&
333	phisfi_val)
334
335	use surfdat_h, only: phisfi, albedodat, &
336	zmea, zstd, zsig, zgam, zthe, z0, &
337	tsurf, co2ice, emis, qsurf
338
339	INTEGER, INTENT(IN) :: ims,ime,jms,jme
340	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,MARS_MODE
341	INTEGER :: i,j,subs,nlast
342	REAL, INTENT(IN ) :: CST_AL, phisfi_val, CST_Z0
343	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
344	M_ALBEDO,M_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
345	REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN ) :: M_GW
346
347	DO j = jps,jpe
348	DO i = ips,ipe
349
350	!-----------------------------------!
351	! 1D subscript for physics "cursor" !
352	!-----------------------------------!
353	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
354
355	!---------------------!
356	! Ground geopotential !
357	!---------------------!
358	phisfi(subs) = phisfi_val
359
360	!---------------!
361	! Ground albedo !
362	!---------------!
363	IF (JULYR .ne. 9999) THEN
364	IF (CST_AL == 0) THEN
365	albedodat(subs)=M_ALBEDO(i,j)
366	ELSE
367	albedodat(subs)=CST_AL
368	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** SET CONSTANT ALBEDO ', albedodat
369	ENDIF
370	ELSE
371	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION albedo: ', CST_AL
372	albedodat(subs)=CST_AL
373	ENDIF
374
375	!-----------------------------------------!
376	! Gravity wave parametrization !
377	! NB: usually 0 in mesoscale applications !
378	!-----------------------------------------!
379	IF (JULYR .ne. 9999) THEN
380	zmea(subs)=M_GW(i,1,j)
381	zstd(subs)=M_GW(i,2,j)
382	zsig(subs)=M_GW(i,3,j)
383	zgam(subs)=M_GW(i,4,j)
384	zthe(subs)=M_GW(i,5,j)
385	ELSE
386	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION GWdrag OFF'
387	zmea(subs)=0.
388	zstd(subs)=0.
389	zsig(subs)=0.
390	zgam(subs)=0.
391	zthe(subs)=0.
392	ENDIF
393
394	!----------------------------!
395	! Variable surface roughness !
396	!----------------------------!
397	IF (JULYR .ne. 9999) THEN
398	IF (CST_Z0 == 0) THEN
399	z0(subs) = M_Z0(i,j)
400	ELSE
401	z0(subs) = CST_Z0
402	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
403	ENDIF
404	ELSE
405	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
406	z0(subs)=CST_Z0
407	ENDIF
408	!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
409	IF (z0(subs) == 0.) THEN
410	IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
411	z0(subs) = 0.01
412	ENDIF
413	!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
414	IF (z0(subs) < 0.) THEN
415	PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
416	PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
417	PRINT *, ' -- or check the constant value set in namelist.input'
418	STOP
419	ENDIF
420
421	!-----------------------------------------------!
422	! Ground temperature, emissivity, CO2 ice cover !
423	!-----------------------------------------------!
424	tsurf(subs) = M_TSURF(i,j)
425	emis(subs) = M_EMISS(i,j)
426	co2ice(subs) = M_CO2ICE(i,j)
427
428	!-------------------!
429	! Tracer at surface !
430	!-------------------!
431	qsurf(subs,:)=0. ! default case
432	SELECT CASE (MARS_MODE)
433	CASE(1)
434	qsurf(subs,2)=M_H2OICE(i,j) !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
435	!! ----- retrocompatible ancienne physique
436	!! ----- [H2O ice is last tracer in qsurf in LMD physics]
437	CASE(2)
438	qsurf(subs,1)=0. !! not coupled with lifting for the moment [non remobilise]
439	!CASE(3)
440	!qsurf(subs,1)=q_prof(1,1) !!! temporaire, a definir
441	!qsurf(subs,2)=q_prof(1,2)
442	CASE(11)
443	qsurf(subs,2)=M_H2OICE(i,j) !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
444	qsurf(subs,3)=0. !! not coupled with lifting for the moment [non remobilise]
445	CASE(12)
446	qsurf(subs,2)=M_H2OICE(i,j) !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
447	qsurf(subs,3)=0. !! not coupled with lifting for the moment [non remobilise]
448	END SELECT
449
450	ENDDO
451	ENDDO
452
453	!!---------------------!!
454	!! OUTPUT FOR CHECKING !!
455	!!---------------------!!
456	nlast = (ipe-ips+1)*(jpe-jps+1)
457	print*,"check: phisfi",phisfi(1),phisfi(nlast)
458	print*,"check: albedodat",albedodat(1),albedodat(nlast)
459	print*,"check: zmea",zmea(1),zmea(nlast)
460	print*,"check: zstd",zstd(1),zstd(nlast)
461	print*,"check: zsig",zsig(1),zsig(nlast)
462	print*,"check: zgam",zgam(1),zgam(nlast)
463	print*,"check: zthe",zthe(1),zthe(nlast)
464	print*,"check: z0",z0(1),z0(nlast)
465	print*,"check: tsurf",tsurf(1),tsurf(nlast)
466	print*,"check: emis",emis(1),emis(nlast)
467	print*,"check: co2ice",co2ice(1),co2ice(nlast)
468	print*,"check: qsurf",qsurf(1,:),qsurf(nlast,:)
469
470	END SUBROUTINE update_inputs_physiq_surf
471
472	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
473	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
474	SUBROUTINE update_inputs_physiq_soil( &
475	ims,ime,jms,jme,&
476	ips,ipe,jps,jpe,&
477	JULYR,nsoil,&
478	M_TI,CST_TI,&
479	M_ISOIL,M_DSOIL,&
480	M_TSOIL,M_TSURF)
481
482	use comsoil_h, only: inertiedat,mlayer,layer,tsoil
483
484	INTEGER, INTENT(IN) :: ims,ime,jms,jme
485	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,nsoil
486	INTEGER :: i,j,subs,nlast,k
487	REAL, INTENT(IN ) :: CST_TI
488	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
489	M_TI, M_TSURF
490	REAL, DIMENSION( ims:ime, nsoil, jms:jme ), INTENT(IN) :: &
491	M_TSOIL, M_ISOIL, M_DSOIL
492	REAL :: inertiedat_val
493	REAL :: lay1,alpha
494
495	DO j = jps,jpe
496	DO i = ips,ipe
497
498	!-----------------------------------!
499	! 1D subscript for physics "cursor" !
500	!-----------------------------------!
501	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
502
503	!-----------------!
504	! Thermal Inertia !
505	!-----------------!
506	IF (JULYR .ne. 9999) THEN
507	IF (CST_TI == 0) THEN
508	inertiedat_val=M_TI(i,j)
509	ELSE
510	inertiedat_val=CST_TI
511	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** SET CONSTANT THERMAL INERTIA ', inertiedat_val
512	ENDIF
513	ELSE
514	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION inertia: ',CST_TI
515	inertiedat_val=CST_TI
516	ENDIF
517	!inertiedat(subs) = inertiedat_val
518	!--pb de dimensions???!!???
519	IF (JULYR .ne. 9999) THEN
520	inertiedat(subs,:)=M_ISOIL(i,:,j) !! verifier que cest bien hires TI en surface
521	mlayer(0:nsoil-1)=M_DSOIL(i,:,j)
522	ELSE
523	IF ( nsoil .lt. 18 ) THEN
524	PRINT ,'* Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
525	STOP
526	ENDIF
527	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
528	do k=1,nsoil
529	inertiedat(subs,k) = inertiedat_val
530	!mlayer(k-1) = sqrt(887.75/3.14)((2.(k-0.5))-1.) inertiedat_val / wvolcapa ! old setting
531	mlayer(k-1) = 2.E-4 * (2.**(k-0.5-1.)) ! new gcm settings
532	enddo
533	ENDIF
534	IF ( (i == ips) .AND. (j == jps) ) &
535	PRINT ,'* Mars ** TI and depth profiles are',inertiedat(subs,:),mlayer(0:nsoil-1)
536
537	!!!!!!!!!!!!!!!!! DONE in soil_setting.F
538	! 1.5 Build layer(); following the same law as mlayer()
539	! Assuming layer distribution follows mid-layer law:
540	! layer(k)=lay1alpha*(k-1)
541	lay1=sqrt(mlayer(0)*mlayer(1))
542	alpha=mlayer(1)/mlayer(0)
543	do k=1,nsoil
544	layer(k)=lay1(alpha*(k-1))
545	enddo
546
547	!------------------------!
548	! Deep soil temperatures !
549	!------------------------!
550	IF (M_TSOIL(i,1,j) .gt. 0. .and. JULYR .ne. 9999) THEN
551	tsoil(subs,:)=M_TSOIL(i,:,j)
552	ELSE
553	IF ( (i == ips) .AND. (j == jps) ) PRINT ,'* Mars ** no tsoil. set it to tsurf.'
554	do k=1,nsoil
555	tsoil(subs,k) = M_TSURF(i,j)
556	enddo
557	ENDIF
558
559	ENDDO
560	ENDDO
561
562	!!---------------------!!
563	!! OUTPUT FOR CHECKING !!
564	!!---------------------!!
565	nlast = (ipe-ips+1)*(jpe-jps+1)
566	print*,"check: inertiedat",inertiedat(1,:),inertiedat(nlast,:)
567	print*,"check: mlayer",mlayer(:)
568	print*,"check: layer",layer(:)
569	print*,"check: tsoil",tsoil(1,:),tsoil(nlast,:)
570
571	END SUBROUTINE update_inputs_physiq_soil
572
573	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
575	SUBROUTINE update_inputs_physiq_turb( &
576	ims,ime,jms,jme,kms,kme,&
577	ips,ipe,jps,jpe,&
578	RESTART,isles,&
579	M_Q2,M_WSTAR)
580
581	use turb_mod, only: q2,wstar,turb_resolved
582
583	INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
584	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
585	INTEGER :: i,j,subs,nlast
586	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: M_WSTAR
587	REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(IN) :: M_Q2
588	LOGICAL, INTENT(IN ) :: RESTART,isles
589
590	!! to know if this is turbulence-resolving run or not
591	turb_resolved = isles
592
593	DO j = jps,jpe
594	DO i = ips,ipe
595
596	!-----------------------------------!
597	! 1D subscript for physics "cursor" !
598	!-----------------------------------!
599	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
600
601	!PBL wind variance
602	IF (.not. restart) THEN
603	q2(subs,:) = 1.E-6
604	wstar(subs)=0.
605	ELSE
606	q2(subs,:)=M_Q2(i,:,j)
607	wstar(subs)=M_WSTAR(i,j)
608	ENDIF
609
610	ENDDO
611	ENDDO
612
613	!!---------------------!!
614	!! OUTPUT FOR CHECKING !!
615	!!---------------------!!
616	nlast = (ipe-ips+1)*(jpe-jps+1)
617	print*,"check: q2",q2(1,:),q2(nlast,:)
618	print*,"check: wstar",wstar(1),wstar(nlast)
619
620	END SUBROUTINE update_inputs_physiq_turb
621
622	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
623	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
624	SUBROUTINE update_inputs_physiq_rad( &
625	ims,ime,jms,jme,&
626	ips,ipe,jps,jpe,&
627	RESTART,&
628	M_FLUXRAD)
629
630	use dimradmars_mod, only: fluxrad
631
632	INTEGER, INTENT(IN) :: ims,ime,jms,jme
633	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
634	INTEGER :: i,j,subs,nlast
635	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: M_FLUXRAD
636	LOGICAL, INTENT(IN ) :: RESTART
637
638	DO j = jps,jpe
639	DO i = ips,ipe
640
641	!-----------------------------------!
642	! 1D subscript for physics "cursor" !
643	!-----------------------------------!
644	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
645
646	! fluxrad_save
647	IF (.not. restart) THEN
648	fluxrad(subs)=0.
649	ELSE
650	fluxrad(subs)=M_FLUXRAD(i,j)
651	ENDIF
652	!! et fluxrad_sky ???!???
653
654	ENDDO
655	ENDDO
656
657	!!---------------------!!
658	!! OUTPUT FOR CHECKING !!
659	!!---------------------!!
660	nlast = (ipe-ips+1)*(jpe-jps+1)
661	print*,"check: fluxrad",fluxrad(1),fluxrad(nlast)
662
663	END SUBROUTINE update_inputs_physiq_rad
664
665	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
666	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
667	SUBROUTINE update_inputs_physiq_slope( &
668	ims,ime,jms,jme,&
669	ips,ipe,jps,jpe,&
670	JULYR,&
671	SLPX,SLPY)
672
673	USE module_model_constants, only: DEGRAD
674	USE slope_mod, ONLY: theta_sl, psi_sl
675
676	INTEGER, INTENT(IN) :: ims,ime,jms,jme
677	INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR
678	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: SLPX,SLPY
679	INTEGER :: i,j,subs,nlast
680
681	DO j = jps,jpe
682	DO i = ips,ipe
683
684	!-----------------------------------!
685	! 1D subscript for physics "cursor" !
686	!-----------------------------------!
687	subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
688
689	!-------------------!
690	! Slope inclination !
691	!-------------------!
692	IF (JULYR .ne. 9999) THEN
693	theta_sl(subs)=atan(sqrt( (1000.SLPX(i,j))2 + (1000.SLPY(i,j))**2 ))
694	theta_sl(subs)=theta_sl(subs)/DEGRAD
695	ELSE
696	IF ( (i == ips) .AND. (j == jps) ) PRINT *,'IDEALIZED SIMULATION: theta_sl = 0'
697	theta_sl(subs)=0.
698	ENDIF
699
700	!-------------------------------------------!
701	! Slope orientation; 0 is north, 90 is east !
702	!-------------------------------------------!
703	IF (JULYR .ne. 9999) THEN
704	psi_sl(subs)=-90.*DEGRAD-atan(SLPY(i,j)/SLPX(i,j))
705	if (SLPX(i,j) .ge. 0.) then
706	psi_sl(subs)=psi_sl(subs)-180.*DEGRAD
707	endif
708	psi_sl(subs)=360.*DEGRAD+psi_sl(subs)
709	psi_sl(subs)=psi_sl(subs)/DEGRAD
710	psi_sl(subs) = MODULO(psi_sl(subs)+180.,360.)
711	ELSE
712	IF ( (i == ips) .AND. (j == jps) ) PRINT *,'IDEALIZED SIMULATION: psi_sl = 0'
713	psi_sl(subs) = 0.
714	ENDIF
715
716	ENDDO
717	ENDDO
718
719	!!---------------------!!
720	!! OUTPUT FOR CHECKING !!
721	!!---------------------!!
722	nlast = (ipe-ips+1)*(jpe-jps+1)
723	print*,"check: theta_sl",theta_sl(1),theta_sl(nlast)
724	print*,"check: psi_sl",psi_sl(1),psi_sl(nlast)
725
726	END SUBROUTINE update_inputs_physiq_slope
727
728	END MODULE update_inputs_physiq_mod
729

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

Download in other formats:

Original Format