Context Navigation

etat0phys_netcdf.f90 @ 5797

Last change on this file since 5797 was 5797, checked in by aborella, 4 months ago
Bugfix for saturation adjustment in cirrus mixing + bugfix for contrails sedimentation + new diagnostics + support for unadjusted contrails
File size: 23.0 KB

Line
1	!
2	! $Id$
3	!
4	MODULE etat0phys
5	!
6	!*******************************************************************************
7	! Purpose: Create physical initial state using atmospheric fields from a
8	! database of atmospheric to initialize the model.
9	!-------------------------------------------------------------------------------
10	! Comments:
11	!
12	! * This module is designed to work for Earth (and with ioipsl)
13	!
14	! * etat0phys_netcdf routine can access to NetCDF data through subroutines:
15	! "start_init_phys" for variables contained in file "ECPHY.nc":
16	! 'ST' : Surface temperature
17	! 'CDSW' : Soil moisture
18	! "start_init_orog" for variables contained in file "Relief.nc":
19	! 'RELIEF' : High resolution orography
20	!
21	! * The land mask and corresponding weights can be:
22	! 1) computed using the ocean mask from the ocean model (to ensure ocean
23	! fractions are the same for atmosphere and ocean) for coupled runs.
24	! File name: "o2a.nc" ; variable name: "OceMask"
25	! 2) computed from topography file "Relief.nc" for forced runs.
26	!
27	! * Allowed values for read_climoz flag are 0, 1 and 2:
28	! 0: do not read an ozone climatology
29	! 1: read a single ozone climatology that will be used day and night
30	! 2: read two ozone climatologies, the average day and night climatology
31	! and the daylight climatology
32	!-------------------------------------------------------------------------------
33	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
34	! I have chosen to use the iml+1 as an argument to this routine and we declare
35	! internaly smaller fields when needed. This needs to be cleared once and for
36	! all in LMDZ. A convention is required.
37	!-------------------------------------------------------------------------------
38
39	USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo
40	USE assert_eq_m, ONLY: assert_eq
41	USE dimphy, ONLY: klon
42	USE conf_dat_m, ONLY: conf_dat2d
43	USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
44	solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, &
45	sollw,sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
46	sig1, ftsol, clwcon, fm_therm, wake_Cstar, pctsrf, entr_therm,radpas, f0,&
47	zmax0,fevap, rnebcon,falb_dir, falb_dif, wake_fip, agesno, detr_therm, pbl_tke, &
48	phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
49	prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
50	ale_bl, ale_bl_trig, alp_bl, &
51	ale_wake, ale_bl_stat, AWAKE_S, &
52	cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, qic_ancien, nic_ancien
53
54	USE comconst_mod, ONLY: pi, dtvr
55	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
56	USE paramet_mod_h
57	USE dimsoil_mod_h, ONLY: nsoilmx
58	USE comgeom2_mod_h
59	USE clesphys_mod_h
60	USE iniprint_mod_h
61	PRIVATE
62	PUBLIC :: etat0phys_netcdf
63
64	REAL, SAVE :: deg2rad
65	REAL, SAVE, ALLOCATABLE :: tsol(:)
66	INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
67	REAL, ALLOCATABLE, SAVE :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
68	CHARACTER(LEN=256), PARAMETER :: oroparam="oro_params.nc"
69	CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
70	CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc", psrfvar="SP"
71	CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW", tsrfvar="ST"
72
73
74	CONTAINS
75
76
77	!-------------------------------------------------------------------------------
78	!
79	SUBROUTINE etat0phys_netcdf(masque, phis)
80	!
81	!-------------------------------------------------------------------------------
82	! Purpose: Creates initial states
83	!-------------------------------------------------------------------------------
84	! Notes: 1) This routine is designed to work for Earth
85	! 2) If masque(:,:)/=-99999., masque and phis are already known.
86	! Otherwise: compute it.
87	!-------------------------------------------------------------------------------
88	USE control_mod
89	USE fonte_neige_mod
90	USE pbl_surface_mod
91	USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
92	USE indice_sol_mod
93	USE conf_phys_m, ONLY: conf_phys
94	USE init_ssrf_m, ONLY: start_init_subsurf
95	USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
96	ratqs_inter_
97	USE alpale_mod
98	USE compbl_mod_h
99	IMPLICIT NONE
100	!-------------------------------------------------------------------------------
101	! Arguments:
102	REAL, INTENT(INOUT) :: masque(:,:) !--- Land mask dim(iip1,jjp1)
103	REAL, INTENT(INOUT) :: phis (:,:) !--- Ground geopotential dim(iip1,jjp1)
104	!-------------------------------------------------------------------------------
105	! Local variables:
106	CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
107	INTEGER :: i, j, l, ji, iml, jml
108	LOGICAL :: read_mask
109	REAL :: phystep, dummy
110	REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
111	REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder
112	!GG
113	REAL, DIMENSION(klon) :: hice, tsic, bilg_cumul
114	!GG
115	REAL, DIMENSION(klon,nbsrf) :: qsurf, snsrf
116	REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
117
118	!--- Arguments for conf_phys
119	LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
120	REAL :: solarlong0, seuil_inversion, fact_cldcon, facttemps
121	LOGICAL :: ok_newmicro
122	INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
123	REAL :: ratqsbas, ratqshaut, tau_ratqs
124	LOGICAL :: ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
125	INTEGER :: flag_aerosol
126	INTEGER :: flag_aerosol_strat
127	INTEGER :: flag_volc_surfstrat
128	LOGICAL :: flag_aer_feedback
129	LOGICAL :: flag_bc_internal_mixture
130	REAL :: bl95_b0, bl95_b1
131	INTEGER :: read_climoz !--- Read ozone climatology
132	REAL :: alp_offset
133	LOGICAL :: filtre_oro=.false.
134
135	deg2rad= pi/180.0
136	iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
137	jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
138
139	! Physics configuration
140	!*******************************************************************************
141	CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
142	callstats, &
143	solarlong0,seuil_inversion, &
144	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
145	iflag_cldcon, &
146	ratqsbas,ratqshaut,tau_ratqs, &
147	!GG iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
148	ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, &
149	aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat, &
150	flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1, &
151	read_climoz, alp_offset)
152	CALL phys_state_var_init(read_climoz)
153
154	!--- Initial atmospheric CO2 conc. from .def file
155	co2_ppm0 = co2_ppm
156
157	! Compute ground geopotential, sub-cells quantities and possibly the mask.
158	!*******************************************************************************
159	read_mask=ANY(masque/=-99999.); masque_tmp=masque
160	CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
161
162	CALL getin('filtre_oro',filtre_oro)
163	IF (filtre_oro) CALL filtreoro(size(phis,1),size(phis,2),phis,masque_tmp,rlatu)
164
165	WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
166	IF(.NOT.read_mask) THEN !--- Keep mask form orography
167	masque=masque_tmp
168	IF(prt_level>=1) THEN
169	WRITE(lunout,*)'BUILT MASK :'
170	WRITE(lunout,fmt) NINT(masque)
171	END IF
172	WHERE( masque(:,:)<EPSFRA) masque(:,:)=0.
173	WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
174	END IF
175	CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
176
177	! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
178	!*******************************************************************************
179	CALL start_init_phys(rlonu, rlatv, phis)
180
181	! Some initializations.
182	!*******************************************************************************
183	sn (:) = 0.0 !--- Snow
184	radsol(:) = 0.0 !--- Net radiation at ground
185	rugmer(:) = 0.001 !--- Ocean rugosity
186	!--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
187	IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
188
189	! Sub-surfaces initialization.
190	!*******************************************************************************
191	CALL start_init_subsurf(read_mask)
192
193	! Write physical initial state
194	!*******************************************************************************
195	WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
196	phystep = dtvr * FLOAT(iphysiq)
197	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
198	WRITE(lunout,*)'phystep =', phystep, radpas
199
200	! Init: ftsol, snsrf, qsurf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
201	!*******************************************************************************
202	DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
203	DO i=1,nbsrf; snsrf(:,i) = sn; END DO
204	falb_dir(:, :, is_ter) = 0.08
205	falb_dir(:, :, is_lic) = 0.6
206	falb_dir(:, :, is_oce) = 0.5
207	falb_dir(:, :, is_sic) = 0.6
208
209	!ym warning missing init for falb_dif => set to 0
210	falb_dif(:,:,:)=0
211
212	u10m(:,:)=0
213	v10m(:,:)=0
214	treedrg(:,:,:)=0
215
216	fevap(:,:) = 0.
217	qsurf = 0.
218	DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
219	rain_fall = 0.
220	snow_fall = 0.
221	solsw = 165.
222	solswfdiff = 1.
223	sollw = -53.
224	!ym warning missing init for sollwdown => set to 0
225	sollwdown = 0.
226	t_ancien = 273.15
227	q_ancien = 0.
228	ql_ancien = 0.
229	qs_ancien = 0.
230	prlw_ancien = 0.
231	prsw_ancien = 0.
232	prw_ancien = 0.
233	agesno = 0.
234
235	u_ancien = 0.
236	v_ancien = 0.
237	wake_delta_pbl_TKE(:,:,:)=0
238	wake_dens(:)=0
239	awake_dens = 0.
240	cv_gen = 0.
241	ale_bl = 0.
242	ale_bl_trig =0.
243	alp_bl=0.
244	ale_wake=0.
245	ale_bl_stat=0.
246
247	cf_ancien = 0.
248	qvc_ancien = 0.
249	qvcon = 0.
250	qccon = 0.
251	cfc_ancien = 0.
252	qtc_ancien = 0.
253	qic_ancien = 0.
254	nic_ancien = 0.
255
256	z0m(:,:)=0 ! ym missing 5th subsurface initialization
257
258	z0m(:,is_oce) = rugmer(:)
259	z0m(:,is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
260	z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
261	z0m(:,is_sic) = 0.001
262	z0h(:,:)=z0m(:,:)
263
264	fder = 0.0
265	clwcon = 0.0
266	rnebcon = 0.0
267	ratqs = 0.0
268	run_off_lic_0 = 0.0
269	rugoro = 0.0
270
271	! Before phyredem calling, surface modules and values to be saved in startphy.nc
272	! are initialized
273	!*******************************************************************************
274	dummy = 1.0
275	pbl_tke(:,:,:) = 1.e-8
276	zmax0(:) = 40.
277	f0(:) = 1.e-5
278	sig1(:,:) = 0.
279	w01(:,:) = 0.
280	wake_deltat(:,:) = 0.
281	wake_deltaq(:,:) = 0.
282	wake_s(:) = 0.
283	wake_cstar(:) = 0.
284	wake_fip(:) = 0.
285	wake_pe = 0.
286	fm_therm = 0.
287	entr_therm = 0.
288	detr_therm = 0.
289	awake_s = 0.
290	!GG
291	hice = 1.0
292	tsic = tsol
293	bilg_cumul = 0.
294	!GG
295
296	CALL fonte_neige_init(run_off_lic_0)
297	!GG CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
298	CALL pbl_surface_init( fder, snsrf, qsurf, tsoil, hice, tsic, bilg_cumul)
299	!GG
300
301	IF (iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1) then
302	delta_tsurf = 0.
303	beta_aridity = 0.
304	end IF
305
306	ratqs_inter_ = 0.002
307	CALL phyredem( "startphy.nc" )
308
309	! WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
310	! WRITE(lunout,*)'entree histclo'
311	CALL histclo()
312
313	END SUBROUTINE etat0phys_netcdf
314	!
315	!-------------------------------------------------------------------------------
316
317
318	!-------------------------------------------------------------------------------
319	!
320	SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
321	!
322	!===============================================================================
323	! Comment:
324	! This routine launch grid_noro, which computes parameters for SSO scheme as
325	! described in LOTT & MILLER (1997) and LOTT(1999).
326	! In case the file oroparam is present and the key read_orop is activated,
327	! grid_noro is bypassed and sub-cell parameters are read from the file.
328	!===============================================================================
329	USE grid_noro_m, ONLY: grid_noro, read_noro
330	USE logic_mod, ONLY: read_orop
331	IMPLICIT NONE
332	!-------------------------------------------------------------------------------
333	! Arguments:
334	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
335	REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
336	!-------------------------------------------------------------------------------
337	! Local variables:
338	CHARACTER(LEN=256) :: modname
339	INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
340	INTEGER :: ierr
341	REAL :: lev(1), date, dt
342	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
343	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var (:,:)
344	REAL, ALLOCATABLE :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
345	REAL, ALLOCATABLE :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
346	!-------------------------------------------------------------------------------
347	modname="start_init_orog"
348	iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
349	jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
350
351	!--- HIGH RESOLUTION OROGRAPHY
352	CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
353
354	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
355	CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
356	lev, ttm_tmp, itau, date, dt, fid)
357	ALLOCATE(relief_hi(iml_rel,jml_rel))
358	CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
359	CALL flinclo(fid)
360
361	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
362	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
363	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
364	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
365
366	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
367	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
368	CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
369	DEALLOCATE(lon_ini,lat_ini)
370
371	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
372	WRITE(lunout,*)
373	WRITE(lunout,)' Compute parameters needed for gravity wave drag code *'
374
375	!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
376	ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
377	ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
378	zsig0(:,:)=0 !ym uninitialized variable
379	zgam0(:,:)=0 !ym uninitialized variable
380	ALLOCATE(zthe0(iml,jml)) !--- Highest slope orientation
381	zthe0(:,:)=0 !ym uninitialized variable
382	ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
383
384	!--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
385	OPEN(UNIT=66,FILE=oroparam,STATUS='OLD',IOSTAT=ierr)
386	IF(ierr==0.AND.read_orop) THEN
387	CLOSE(UNIT=66)
388	CALL read_noro(lon_in,lat_in,oroparam, &
389	phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
390	ELSE
391	!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
392	CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in, &
393	phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
394	END IF
395	phis = phis * 9.81
396	phis(iml,:) = phis(1,:)
397	DEALLOCATE(relief_hi,lon_rad,lat_rad)
398
399	!--- PUT QUANTITIES TO PHYSICAL GRID
400	CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
401	CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
402	CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
403	CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
404	CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
405	CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
406	CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
407
408
409	END SUBROUTINE start_init_orog
410	!
411	!-------------------------------------------------------------------------------
412
413
414	!-------------------------------------------------------------------------------
415	!
416	SUBROUTINE start_init_phys(lon_in,lat_in,phis)
417	!
418	!===============================================================================
419	! Purpose: Compute tsol and qsol, knowing phis.
420	!===============================================================================
421	IMPLICIT NONE
422	!-------------------------------------------------------------------------------
423	! Arguments:
424	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml2)
425	REAL, INTENT(IN) :: phis(:,:) ! dim (iml,jml)
426	!-------------------------------------------------------------------------------
427	! Local variables:
428	CHARACTER(LEN=256) :: modname
429	REAL :: date, dt
430	INTEGER :: iml, jml, jml2, itau(1)
431	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:)
432	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
433	REAL, ALLOCATABLE :: ts(:,:), qs(:,:)
434	!-------------------------------------------------------------------------------
435	modname="start_init_phys"
436	iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
437	jml=SIZE(phis,2); jml2=SIZE(lat_in)
438
439	WRITE(lunout,*)'Opening the surface analysis'
440	CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
441	WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys
442
443	ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
444	ALLOCATE(levphys_ini(llm_phys))
445	CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys, &
446	lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
447
448	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
449	ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
450	lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
451	lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
452
453	ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
454	CALL get_var_phys(tsrfvar,ts) !--- SURFACE TEMPERATURE
455	CALL get_var_phys(qsolvar,qs) !--- SOIL MOISTURE
456	CALL flinclo(fid_phys)
457	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
458
459	!--- TSOL AND QSOL ON PHYSICAL GRID
460	ALLOCATE(tsol(klon))
461	CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
462	CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
463	DEALLOCATE(ts,qs)
464
465	CONTAINS
466
467	!-------------------------------------------------------------------------------
468	!
469	SUBROUTINE get_var_phys(title,field)
470	!
471	!-------------------------------------------------------------------------------
472	IMPLICIT NONE
473	!-------------------------------------------------------------------------------
474	! Arguments:
475	CHARACTER(LEN=*), INTENT(IN) :: title
476	REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
477	!-------------------------------------------------------------------------------
478	! Local variables:
479	INTEGER :: tllm
480	!-------------------------------------------------------------------------------
481	SELECT CASE(title)
482	CASE(psrfvar); tllm=0
483	CASE(tsrfvar,qsolvar); tllm=llm_phys
484	END SELECT
485	IF(ALLOCATED(field)) RETURN
486	ALLOCATE(field(iml,jml)); field(:,:)=0.
487	CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
488	CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
489	CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, &
490	lon_in, lat_in, field)
491
492	END SUBROUTINE get_var_phys
493	!
494	!-------------------------------------------------------------------------------
495	!
496	END SUBROUTINE start_init_phys
497	!
498	!-------------------------------------------------------------------------------
499
500
501	!-------------------------------------------------------------------------------
502	!
503	SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
504	!
505	!-------------------------------------------------------------------------------
506	USE inter_barxy_m, ONLY: inter_barxy
507	IMPLICIT NONE
508	!-------------------------------------------------------------------------------
509	! Arguments:
510	CHARACTER(LEN=*), INTENT(IN) :: nam
511	LOGICAL, INTENT(IN) :: ibeg
512	REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj)
513	REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj)
514	REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
515	REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1)
516	!-------------------------------------------------------------------------------
517	! Local variables:
518	CHARACTER(LEN=256) :: modname
519	INTEGER :: ii, jj, i1, j1, j2
520	REAL, ALLOCATABLE :: vtmp(:,:)
521	!-------------------------------------------------------------------------------
522	modname="interp_startvar"
523	ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
524	jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
525	i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
526	j1=SIZE(varo,2); j2=SIZE(lat2)
527	ALLOCATE(vtmp(i1-1,j1))
528	IF(ibeg.AND.prt_level>1) THEN
529	WRITE(lunout,*)"--------------------------------------------------------"
530	WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
531	WRITE(lunout,*)"--------------------------------------------------------"
532	END IF
533	CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
534	CALL gr_int_dyn(vtmp, varo, i1-1, j1)
535
536	END SUBROUTINE interp_startvar
537	!
538	!-------------------------------------------------------------------------------
539	!
540	!*******************************************************************************
541
542	SUBROUTINE filtreoro(imp1,jmp1,phis,masque,rlatu)
543
544	IMPLICIT NONE
545
546	INTEGER imp1,jmp1
547	REAL, DIMENSION(imp1,jmp1) :: phis,masque
548	REAL, DIMENSION(jmp1) :: rlatu
549	REAL, DIMENSION(imp1) :: wwf
550	REAL, DIMENSION(imp1,jmp1) :: phiso
551	INTEGER :: ifiltre,ifi,ii,i,j
552	REAL :: coslat0,ssz
553
554	coslat0=0.5
555	phiso=phis
556	do j=2,jmp1-1
557	print*,'avant if ',cos(rlatu(j)),coslat0
558	if (cos(rlatu(j))<coslat0) then
559	! nb de pts affectes par le filtrage de part et d'autre du pt
560	ifiltre=(coslat0/cos(rlatu(j))-1.)/2.
561	wwf=0.
562	do i=1,ifiltre
563	wwf(i)=1.
564	enddo
565	wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre
566	do i=1,imp1-1
567	if (masque(i,j)>0.9) then
568	ssz=phis(i,j)
569	do ifi=1,ifiltre+1
570	ii=i+ifi
571	if (ii>imp1-1) ii=ii-imp1+1
572	ssz=ssz+wwf(ifi)*phis(ii,j)
573	ii=i-ifi
574	if (ii<1) ii=ii+imp1-1
575	ssz=ssz+wwf(ifi)*phis(ii,j)
576	enddo
577	phis(i,j)=ssz*cos(rlatu(j))/coslat0
578	endif
579	enddo
580	print,'j=',j,coslat0/cos(rlatu(j)), (1.+2.sum(wwf))*cos(rlatu(j))/coslat0
581	endif
582	enddo
583	call dump2d(imp1,jmp1,phis,'phis ')
584	call dump2d(imp1,jmp1,masque,'masque ')
585	call dump2d(imp1,jmp1,phis-phiso,'dphis ')
586
587	END SUBROUTINE filtreoro
588
589
590	END MODULE etat0phys

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

Download in other formats:

Original Format