Context Navigation

regr_pr_time_av_m.F90 @ 5305

Last change on this file since 5305 was 3598, checked in by jghattas, 5 years ago
Correction sur le commit [3596]. Desole.
File size: 34.7 KB

Line
1	! $Id$
2	MODULE regr_pr_time_av_m
3
4	USE write_field_phy
5	USE mod_phys_lmdz_transfert_para, ONLY: bcast
6	USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
7	USE print_control_mod, ONLY: prt_level
8	IMPLICIT NONE
9
10	!-------------------------------------------------------------------------------
11	! Purpose: Regrid pressure with an averaging method. Operations done:
12	! * on the root process: read a NetCDF 3D or 4D field at a given day.
13	! * pack the fields to the LMDZ "horizontal "physics" grid and scatter
14	! to all threads of all processes;
15	! * in all the threads of all the processes, regrid the fields in pressure
16	! to the LMDZ vertical grid.
17	! * the forcing fields are stretched if the following arguments are present:
18	! - "lat_in": input file latitudes.
19	! - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
20	! so that the tropopause is shifted to the position of the one of LMDZ.
21	! Note that the ozone quantity conservation is not ensured.
22	!-------------------------------------------------------------------------------
23	! Initial routine: regr_pr_av_m module (L. Guez).
24	! Present version: David Cugnet ; corresponding additions:
25	! - time interpolation
26	! - 3D compliant
27	! - vertical stretching of the field to allow tropopause matching between
28	! input field and current lmdz field.
29	!-------------------------------------------------------------------------------
30	! Remarks:
31	! * 3D fields are zonal means, with dimensions (latitude, pressure, julian day)
32	! * 4D fields have the dimensions: (longitude, latitude, pressure, julian day)
33	! * All the fields are already on the horizontal grid (rlatu) or (rlatu,rlonv),
34	! except that the latitudes are in ascending order in the input file.
35	! * We assume that all the input fields have the same coordinates.
36	! * The target vertical LMDZ grid is the grid of layer centers.
37	! * Regridding in pressure can be:
38	! - Ploc=='I': pressures provided at Interfaces => conservative 2nd order
39	! OK for ozone coefficients regridding in Cariolle routines.
40	! - Ploc=='C': pressures provides at cells Centers => classical linear
41	! OK for ozone vertical regridding, especially when torpopause
42	! adjustment is activated, to avoid "strairs shape effect" on profiles.
43	! * All the fields are regridded as a single multi-dimensional array, so it
44	! saves CPU time to call this procedure once for several NetCDF variables
45	! rather than several times, each time for a single NetCDF variable.
46	! * The input file pressure at tropopause can be (in decreasing priority):
47	! 1) read from the file if "tropopause_air_pressure" field is available.
48	! 2) computed using "tro3" and "tro3_at_tropopause' (if available).
49	! 3) computed using "tro3" and a fixed threshold otherwise, constant or
50	! determined using an empirical three parameters law:
51	! o3t(ppbV)=co1+co2SIN(PI(month-2)/6)*TANH(lat_deg/co3)
52	! => co1 and co2 are in ppbV, and co3 in degrees.
53	! => co3 allow a smooth transition between north and south hemispheres.
54	! * If available, the field "ps" (input file pressure at surface) is used.
55	! If not, the current LMDZ ground pressure is taken instead.
56	! * Fields with suffix "m"/"p" are at the closest records earlier/later than
57	! the mid-julian day "julien", on the global "dynamics" horizontal grid.
58	! * Fields(i,j,k,l) are at longitude-latitude-name "rlonv(i)-rlatu(j)-nam(l)",
59	! pressure level/interval (Ploc=="C"/"I") "pcen_in(k)"/"pcen_in(k:k+1)]".
60	! * In the 2D file case, the values are the same for all longitudes.
61	! * The tropopause correction works like this: the input fields (file) are
62	! interpolated on output (LMDZ) pressure field, which is streched using a power
63	! law in a limited zone made of 2 layers:
64	! 1) between lower bound (lower than lowest tropopause) and LMDZ tropopause
65	! 2) between LMDZ tropopause and upper bound (higher thzn highest tropopause)
66	! The stretching function has the following form:
67	! Sigma_str = Sigma^(1+alpha*phi(Sigma)), where:
68	! * alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
69	! This value shifts the file tropopause to the height of the one of LMDZ.
70	! * phi is quasi-linear (sections of 1/log function) in the adjacent intervals:
71	! - from 0 to 1 in [Sig_top,SigT_ou]
72	! - from 1 to 0 in [SigT_ou,Sig_bot]
73	! This quasi-triangular localization function ponderates alpha-law from one near
74	! the tropopause to zero each side apart.
75	!
76	! * The following fields are on the global "dynamics" grid, as read from files:
77	REAL, SAVE, ALLOCATABLE :: v1 (:,:,:,:) !--- Current time ozone fields
78	! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
79	REAL, SAVE, ALLOCATABLE :: v1m(:,:,:,:) !--- Previous time ozone fields
80	REAL, SAVE, ALLOCATABLE :: v1p(:,:,:,:) !--- Next time ozone fields
81	REAL, SAVE, ALLOCATABLE :: pgm(:,:), pgp(:,:) !--- Ground pressure
82	REAL, SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
83	REAL, SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
84	INTEGER, SAVE :: ntim_in !--- Records nb in input file
85	INTEGER, SAVE :: itrp0 !--- idx above chem tropop.
86	INTEGER, SAVE :: irec !--- Current time index
87	! * for daily input files: current julian day number
88	! * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
89	LOGICAL, SAVE :: linterp !--- Interpolation in time
90	LOGICAL, SAVE :: lPrSfile !--- Surface pressure flag
91	LOGICAL, SAVE :: lPrTfile !--- Tropopause pressure flag
92	LOGICAL, SAVE :: lO3Tfile !--- Tropopause ozone flag
93	LOGICAL, SAVE :: lfirst=.TRUE. !--- First call flag
94	!$OMP THREADPRIVATE(lfirst)
95	REAL, PARAMETER :: pTropUp=9.E+3 !--- Value < tropopause pressure (Pa)
96	REAL, PARAMETER :: gamm =0.4 !--- Max. stretched layer sigma ratio
97	REAL, PARAMETER :: rho =1.4 !--- Max tropopauses sigma ratio
98	REAL, PARAMETER :: o3t0 =1.E-7 !--- Nominal O3 vmr at tropopause
99	LOGICAL, PARAMETER :: lO3Tpara=.FALSE. !--- Parametrized O3 vmr at tropopause
100	LOGICAL, PARAMETER :: ldebug=.FALSE.!--- Force writefield_phy multiple outputs
101	REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum
102	REAL, PARAMETER :: ChemPTrMax=4.E+4 ! chemical tropopause pressure (Pa).
103	REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum
104	REAL, PARAMETER :: DynPTrMax =4.E+4 ! dynamical tropopause pressure (Pa).
105
106	CONTAINS
107
108	!-------------------------------------------------------------------------------
109	!
110	SUBROUTINE regr_pr_time_av(fID, nam, julien, Ploc, Pre_in, Pre_ou, v3, Pgnd_ou,&
111	time_in, lon_in, lat_in, Ptrp_ou)
112	!
113	!-------------------------------------------------------------------------------
114	USE dimphy, ONLY: klon
115	USE netcdf95, ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
116	NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
117	USE netcdf, ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
118	USE assert_m, ONLY: assert
119	USE assert_eq_m, ONLY: assert_eq
120	!! USE comvert_mod, ONLY: scaleheight
121	USE interpolation, ONLY: locate
122	USE regr_conserv_m, ONLY: regr_conserv
123	USE regr_lint_m, ONLY: regr_lint
124	USE slopes_m, ONLY: slopes
125	USE mod_phys_lmdz_para, ONLY: is_mpi_root,is_master
126	USE mod_grid_phy_lmdz, ONLY: nlon=>nbp_lon, nlat=>nbp_lat, nlev_ou=>nbp_lev, klon_glo, grid_type, unstructured
127	USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter, gather
128	USE phys_cal_mod, ONLY: calend, year_len, days_elapsed, jH_cur
129	USE geometry_mod, ONLY: ind_cell_glo
130	!-------------------------------------------------------------------------------
131	! Arguments:
132	INTEGER, INTENT(IN) :: fID !--- NetCDF file ID
133	CHARACTER(LEN=13), INTENT(IN) :: nam(:) !--- NetCDF variables names
134	REAL, INTENT(IN) :: julien !--- Days since Jan 1st
135	CHARACTER(LEN=1), INTENT(IN) :: Ploc !--- Pressures locations
136	!--- File/LMDZ (resp. decreasing & increasing order) pressure, Pa
137	! At cells centers or interfaces depending on "Ploc" keyword (C/I)
138	REAL, INTENT(IN) :: Pre_in(:) !--- in: file (nlev_in[+1])
139	REAL, INTENT(IN) :: Pre_ou(:,:) !--- out: LMDZ (klon,nlev_ou[+1])
140	REAL, INTENT(OUT) :: v3(:,:,:) !--- Regr. fld (klon,nlev_ou,n_var)
141	REAL, INTENT(IN), OPTIONAL :: Pgnd_ou(:) !--- LMDZ ground pressure (klon)
142	REAL, INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
143	! since Jan 1 of current year
144	REAL, INTENT(IN), OPTIONAL :: lon_in(:) !--- File longitudes vector (klon)
145	REAL, INTENT(IN), OPTIONAL :: lat_in(:) !--- File latitudes vector (klon)
146	REAL, INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres (klon)
147	!-------------------------------------------------------------------------------
148	! Local variables:
149	include "clesphys.h"
150	include "YOMCST.h"
151	CHARACTER(LEN=80) :: sub
152	CHARACTER(LEN=320) :: str
153	INTEGER :: vID, ncerr, n_var, ibot, iout, nn
154	INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
155	LOGICAL :: lAdjTro !--- Need to adjust tropopause
156	REAL :: y_frac !--- Elapsed year fraction
157	REAL :: alpha, beta, al !--- For stretching/interpolation
158	REAL :: SigT_in, SigT_ou !--- Input and output tropopauses
159	REAL :: Sig_bot, Sig_top !--- Bounds of quasi-hat function
160	REAL :: Sig_bo0, Sig_to0 !--- Lower/upper tropopauses
161	REAL :: Sig_in (SIZE(Pre_in)) !--- Input field sigma levels
162	REAL :: Sig_ou (SIZE(Pre_ou,2)) !--- Output LMDZ sigma levels
163	REAL :: phi (SIZE(Pre_ou,2)) !--- Stretching exponent anomaly
164	REAL :: Pstr_ou(SIZE(Pre_ou,2)) !--- Stretched pressure levels
165	REAL :: Pres_ou(SIZE(Pre_ou,2)) !--- Pre_ou(i,:), reversed order
166	REAL, DIMENSION(nlon, nlat) :: pg1, pt1, ot1
167	REAL, DIMENSION(klon) :: Pgnd_in, Ptrp_in, Otrp_in
168	REAL, DIMENSION(klon) :: Ptrop_ou, Pgrnd_ou
169	! * The following fields are scattered to the partial "physics" horizontal grid.
170	REAL, POINTER :: v2(:,:,:) !--- Current time ozone fields
171	! In the 2D file case, the values are the same for all longitudes.
172	! "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)" and name "nam(l)"
173	! Both are: * if Ploc=='I' in pressure interval "press_in_edg(k:k+1)"
174	! * if Ploc=='C' at pressure "press_in_cen(k)"
175	REAL, TARGET :: &
176	v2i(klon,SIZE(Pre_in)-1,SIZE(nam)), & !--- v2 in Ploc=='I' case
177	v2c(klon,SIZE(Pre_in) ,SIZE(nam)) !--- v2 in Ploc=='C' case
178	INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
179	LOGICAL :: ll
180	!--- For debug
181	REAL, DIMENSION(klon) :: Ptrop_in, Ptrop_ef
182	REAL, DIMENSION(klon) :: dzStrain, dzStrain0
183	REAL, DIMENSION(klon,SIZE(Pre_ou,2)) :: Pstrn_ou, phii
184	!-------------------------------------------------------------------------------
185	sub="regr_pr_time_av"
186	nlev_in=SIZE(Pre_in); IF(Ploc=='I') nlev_in=nlev_in-1
187	IF(Ploc=='I') THEN; v2 => v2i; ELSE; v2 => v2c; END IF
188	CALL assert(SIZE(Pre_ou,1)==klon,TRIM(sub)//" Pre_ou klon")
189	CALL assert(SIZE(v3,1)==klon, TRIM(sub)//" v3 klon")
190	CALL assert(SIZE(v3,2)==nlev_ou, TRIM(sub)//" v3 nlev_ou")
191	IF(Ploc=='I') CALL assert(SIZE(Pre_ou,2)==nlev_ou+1,TRIM(sub)//" Pre_ou nlev_ou+1")
192	IF(Ploc=='C') CALL assert(SIZE(Pre_ou,2)==nlev_ou ,TRIM(sub)//" Pre_ou nlev_ou")
193	n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
194	IF(PRESENT(Pgnd_ou)) CALL assert(SIZE(Pgnd_ou)==klon,TRIM(sub)//" Pgnd_ou klon")
195	IF(PRESENT(lon_in)) CALL assert(SIZE(lon_in )==klon,TRIM(sub)//" lon_in klon")
196	IF(PRESENT(lat_in)) CALL assert(SIZE(lat_in )==klon,TRIM(sub)//" lat_in klon")
197	IF(PRESENT(Ptrp_ou)) CALL assert(SIZE(Ptrp_ou)==klon,TRIM(sub)//" Ptrp_ou klon")
198	lAdjTro=PRESENT(Ptrp_ou)
199	IF(lAdjTro) THEN
200	IF(.NOT.PRESENT(lat_in)) &
201	CALL abort_physic(sub, 'Missing lat_in (required if adjust_tropopause=T)', 1)
202	IF(.NOT.PRESENT(Pgnd_ou).AND.Ploc=='C') &
203	CALL abort_physic(sub, 'Missing ground Pr(required if adjust_tropopause=T)', 1)
204	IF(PRESENT(Pgnd_ou)) THEN; Pgrnd_ou=Pgnd_ou; ELSE; Pgrnd_ou=Pre_ou(:,1); END IF
205	END IF
206
207	!$OMP MASTER
208	IF(is_mpi_root) THEN
209
210	!=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
211	IF(lfirst) THEN
212	lPrSfile=lAdjTro.AND.NF90_INQ_VARID(fID,"ps" ,vID)==NF90_NOERR
213	lPrTfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
214	lO3Tfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tro3_at_tropopause" ,vID)==NF90_NOERR
215	CALL NF95_INQ_DIMID(fID,"time",vID)
216	CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
217	linterp=PRESENT(time_in).AND.ntim_in==14
218	ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
219	IF(linterp) THEN
220	ALLOCATE(v1m(nlon,nlat,nlev_in,n_var),v1p(nlon,nlat,nlev_in,n_var))
221	IF(lPrSfile) ALLOCATE(pgm(nlon,nlat),pgp(nlon,nlat))
222	IF(lPrTfile) ALLOCATE(ptm(nlon,nlat),ptp(nlon,nlat))
223	IF(lO3Tfile) ALLOCATE(otm(nlon,nlat),otp(nlon,nlat))
224	END IF
225	!--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
226	IF(lAdjTro) itrp0=locate(Pre_in,pTropUp)
227	CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.' ,sub)
228	CALL msg(lPrSfile,'Using GROUND PRESSURE from input O3 forcing file.',sub)
229	CALL msg(lAdjTro ,'o3 forcing file tropopause location uses:' ,sub)
230	IF(lPrTfile) THEN; str=' INPUT FILE PRESSURE'
231	ELSE IF(lO3Tfile) THEN; str=' INPUT FILE O3 CONCENTRATION'
232	ELSE IF(lO3Tpara) THEN; str=' PARAMETRIZED O3 concentration'
233	ELSE; str=' CONSTANT O3 concentration'; END IF
234	CALL msg(lAdjTro,TRIM(str)//' at tropopause')
235	END IF
236
237	!=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
238	CALL update_fields()
239
240	!=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
241	IF(linterp) THEN
242	WRITE(str,'(a,f12.8,2(a,f5.1))')'Interpolating O3 at julian day ',julien,&
243	' from fields at times ',time_in(irec),' and ', time_in(irec+1)
244	CALL msg(.TRUE.,str,sub)
245	al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
246	v1=alv1m+(1.-al)v1p
247	IF(lPrSfile) pg1=alpgm+(1.-al)pgp
248	IF(lPrTfile) pt1=alptm+(1.-al)ptp
249	IF(lO3Tfile) ot1=alotm+(1.-al)otp
250	END IF
251	ELSE
252	IF (lfirst) ALLOCATE(v1(0,0,0,0))
253	END IF
254	!$OMP END MASTER
255	!$OMP BARRIER
256	IF(lfirst) THEN
257	lfirst=.FALSE.; CALL bcast(lfirst)
258	IF(lAdjTro) CALL bcast(itrp0)
259	CALL bcast(lPrSfile); CALL bcast(lPrTfile)
260	CALL bcast(lO3Tfile); CALL bcast(linterp)
261	END IF
262
263	IF (is_master) THEN
264	ALLOCATE(ind_cell_glo_glo(klon_glo))
265	ELSE
266	ALLOCATE(ind_cell_glo_glo(0))
267	ENDIF
268	CALL gather(ind_cell_glo,ind_cell_glo_glo)
269	IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
270
271	CALL scatter2d(v1,v2)
272
273	!--- No "ps" in input file => assumed to be equal to current LMDZ ground press
274	IF(lPrSfile) THEN
275	IF (is_master .AND. grid_type==unstructured) pg1(:,:)=pg1(:,ind_cell_glo_glo(:))
276	CALL scatter2d(pg1,Pgnd_in)
277	ELSE
278	Pgnd_in=Pre_ou(:,1)
279	END IF
280
281	IF(lPrTfile) THEN
282	IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
283	CALL scatter2d(pt1,Ptrp_in)
284	ENDIF
285
286	IF(lO3Tfile) THEN
287	IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
288	CALL scatter2d(ot1,Otrp_in)
289	ENDIF
290	!--- No ground pressure in input file => choose it to be the one of LMDZ
291	IF(lAdjTro.AND..NOT.lPrSfile) Pgnd_in(:)=Pgrnd_ou(:)
292
293	!--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
294	IF(.NOT.lAdjTro) THEN
295	DO i=1,klon
296	Pres_ou=Pre_ou(i,SIZE(Pre_ou,2):1:-1) !--- pplay & paprs are decreasing
297	IF(Ploc=='C') CALL regr_lint (1,v2(i,:,:), LOG(Pre_in(:)), &
298	LOG(Pres_ou(:)), v3(i,nlev_ou:1:-1,:))
299	IF(Ploc=='I') CALL regr_conserv(1,v2(i,:,:), Pre_in(:) , &
300	Pres_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
301	END DO
302	!-------------------------------------------------------------------------------
303	ELSE !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
304	!-------------------------------------------------------------------------------
305	y_frac=(REAL(days_elapsed)+jH_cur)/year_len
306
307	!--- OUTPUT SIGMA LEVELS
308	DO i=1,klon
309
310	!--- INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN
311	Pres_ou = Pre_ou(i,SIZE(Pre_ou,2):1:-1)!--- pplay & paprs are decreasing
312	Sig_in(:) = Pre_in (:)/Pgnd_in(i) !--- increasing values
313	Sig_ou(:) = Pres_ou(:)/Pgnd_ou(i) !--- increasing values
314
315	!--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
316	! to keep tropopause pressure realistic ; high values are usually due to
317	! ozone hole fooling the crude chemical tropopause detection algorithm.
318	SigT_in = get_SigTrop(i,itrp)
319	SigT_in=MIN(SigT_in,ChemPTrMax/Pgnd_in(i)) !--- too low value filtered
320	SigT_in=MAX(SigT_in,ChemPTrMin/Pgnd_ou(i)) !--- too high value filtered
321
322	!--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
323	! dynamical tropopause (especially in filaments) are conterbalanced with
324	! a filter ensuring it stays within a certain distance around input (file)
325	! tropopause, hence avoiding avoid a too thick stretched region ; a final
326	! extra-safety filter keeps the tropopause pressure value realistic.
327	SigT_ou = Ptrp_ou(i)/Pgnd_ou(i)
328	IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low value w/r input
329	IF(SigT_ou>SigT_inrho) SigT_ou=SigT_inrho !--- too high value w/r input
330	SigT_ou=MIN(SigT_ou,DynPTrMax/Pgnd_ou(i)) !--- too low value filtered
331	SigT_ou=MAX(SigT_ou,DynPTrMin/Pgnd_ou(i)) !--- too high value filtered
332	Ptrop_ou(i)=SigT_ou*Pgnd_ou(i)
333	iout = locate(Sig_ou(:),SigT_ou)
334
335	!--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
336	alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
337
338	!--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS
339	Sig_bo0 = MAX(SigT_in,SigT_ou) !--- lowest tropopause
340	Sig_to0 = MIN(SigT_in,SigT_ou) !--- highest tropopause
341	beta = (Sig_bo0/Sig_to0)**gamm !--- stretching exponent
342	Sig_bot = MIN(Sig_bo0beta,0.1(9.+Sig_bo0)) !--- must be <1
343	ibot = locate(Sig_ou(:),Sig_bot) !--- layer index
344	IF(ibot-iout<2) THEN !--- at least one layer thick
345	ibot=MIN(iout+2,nlev_ou); Sig_bot=Sig_ou(ibot)
346	END IF
347	Sig_top = Sig_to0/beta !--- upper bound
348	itop = locate(Sig_ou(:),Sig_top) !--- layer index
349	IF(iout-itop<2) THEN !--- at least one layer thick
350	itop=MAX(iout-2,1); Sig_top=Sig_ou(itop)
351	END IF
352
353	!--- STRETCHING POWER LAW LOCALIZATION FUNCTION:
354	! 0 in [0,Sig_top] 0->1 in [Sig_top,SigT_ou]
355	! 0 in [Sig_bot,1] 1->0 in [SigT_ou, Sig_bot]
356	phi(:)=0.
357	phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&
358	*LOG(SigT_ou)/LOG(SigT_ou/Sig_top)
359	phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&
360	*LOG(SigT_ou)/LOG(SigT_ou/Sig_bot)
361
362	!--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
363	Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)*(alphaphi(:))
364
365	!--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
366	IF(Ploc=='C') CALL regr_lint (1, v2(i,:,:), LOG(Pre_in(:)), &
367	LOG(Pstr_ou(:)), v3(i,nlev_ou:1:-1,:))
368	IF(Ploc=='I') CALL regr_conserv(1, v2(i,:,:), Pre_in(:) , &
369	Pstr_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
370
371	!--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
372	i0=nlev_ou-locate(Pres_ou(:),Ptrop_ou(i))+1
373	ll=check_ozone(v3(i, 1:i0-1 ,1),lon_in(i),lat_in(i),1 ,'troposphere', &
374	5.E-9,3.0E-7)
375	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
376	ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
377	5.E-8,1.5E-5)
378	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
379
380	IF(ldebug) THEN
381	dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)
382	dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)
383	Ptrop_in (i) = SigT_in*Pgnd_in(i)
384	Pstrn_ou(i,:)= Pstr_ou
385	phii(i,:) = phi(:)
386	Ptrop_ef(i) = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp), &
387	Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)
388	END IF
389	END DO
390	END IF
391	IF(ldebug.AND.lAdjTro) THEN
392	CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres
393	CALL writefield_phy('dzStrain' ,dzStrain ,1) !--- Strained thickness
394	CALL writefield_phy('dzStrain0',dzStrain0,1) !--- Tropopauses distance
395	CALL writefield_phy('phi',phii,nlev_ou) !--- Localization function
396	!--- Tropopauses pressures:
397	CALL writefield_phy('PreTr_in',Ptrop_in,1) !--- Input and effective
398	CALL writefield_phy('PreTr_ou',Ptrop_ou,1) !--- LMDz dyn tropopause
399	CALL writefield_phy('PreTr_ef',Ptrop_ef,1) !--- Effective chem tropop
400	END IF
401	IF(ldebug) THEN
402	CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
403	CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field
404	CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure
405	END IF
406
407	CONTAINS
408
409
410	!-------------------------------------------------------------------------------
411	!
412	SUBROUTINE update_fields()
413	!
414	!-------------------------------------------------------------------------------
415	IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION
416	CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
417	irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number
418	!--- MIN -> Security in the unlikely case of roundup errors.
419	CALL get_3Dfields(v1) !--- Read ozone field(s)
420	IF(lAdjTro) THEN !--- Additional files for fields strain
421	IF(lPrSfile) CALL get_2Dfield(pg1,"ps")
422	IF(lPrTfile) CALL get_2Dfield(pt1,"tropopause_air_pressure")
423	IF(lO3Tfile) CALL get_2Dfield(ot1,"tro3_at_tropopause")
424	END IF
425	ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS
426	IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
427	IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
428	CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
429	IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS
430	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&
431	//'reading record ',irec,' (time ',time_in(irec),')'
432	CALL msg(.TRUE.,str,sub)
433	CALL get_3Dfields(v1m) !--- Read ozone field(s)
434	IF(lAdjTro) THEN !--- Additional files for fields strain
435	IF(lPrSfile) CALL get_2Dfield(pgm,"ps")
436	IF(lPrTfile) CALL get_2Dfield(ptm,"tropopause_air_pressure")
437	IF(lO3Tfile) CALL get_2Dfield(otm,"tro3_at_tropopause")
438	END IF
439	ELSE !=== SHIFT FIELDS
440	irec=irec+1
441	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&
442	//' current next one (',irec,', time ',time_in(irec),')'
443	CALL msg(.TRUE.,str,sub)
444	v1m=v1p !--- Ozone fields
445	IF(lAdjTro) THEN !--- Additional files for fields strain
446	IF(lPrSfile) pgm=pgp !--- Surface pressure
447	IF(lPrTfile) ptm=ptp !--- Tropopause pressure
448	IF(lO3Tfile) otm=otp !--- Tropopause ozone
449	END IF
450	END IF
451	irec=irec+1
452	WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'&
453	,irec,' (time ',time_in(irec),')'
454	CALL msg(.TRUE.,str,sub)
455	CALL get_3Dfields(v1p) !--- Read ozone field(s)
456	IF(lAdjTro) THEN !--- Additional files for fields strain
457	IF(lPrSfile) CALL get_2Dfield(pgp,"ps")
458	IF(lPrTfile) CALL get_2Dfield(ptp,"tropopause_air_pressure")
459	IF(lO3Tfile) CALL get_2Dfield(otp,"tro3_at_tropopause")
460	END IF
461	irec=irec-1
462	END IF
463
464	END SUBROUTINE update_fields
465	!
466	!-------------------------------------------------------------------------------
467
468
469	!-------------------------------------------------------------------------------
470	!
471	SUBROUTINE get_2Dfield(v,var)
472	!
473	!-------------------------------------------------------------------------------
474	! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
475	! from the input file.
476	! Remark: In case the field is zonal, it is duplicated along first dimension.
477	!-------------------------------------------------------------------------------
478	! Arguments:
479	REAL, INTENT(INOUT) :: v(:,:)
480	CHARACTER(LEN=*), INTENT(IN) :: var
481	!-------------------------------------------------------------------------------
482	CALL NF95_INQ_VARID(fID, TRIM(var), vID)
483	CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
484	IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[ 1,irec])
485	IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
486	CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
487
488	!--- Flip latitudes: ascending in input file, descending in "rlatu".
489	IF(n_dim==3) THEN
490	v(1,:) = v(1,nlat:1:-1)
491	v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1) !--- Duplication
492	ELSE
493	v(:,:) = v(:,nlat:1:-1)
494	END IF
495
496	END SUBROUTINE get_2Dfield
497	!
498	!-------------------------------------------------------------------------------
499
500
501	!-------------------------------------------------------------------------------
502	!
503	SUBROUTINE get_3Dfields(v)
504	!
505	!-------------------------------------------------------------------------------
506	! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
507	! from the input file. Fields are stacked on fourth dimension.
508	! Remark: In case the field is zonal, it is duplicated along first dimension.
509	!-------------------------------------------------------------------------------
510	! Arguments:
511	REAL, INTENT(INOUT) :: v(:,:,:,:)
512	!-------------------------------------------------------------------------------
513	DO i=1,SIZE(nam)
514	CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
515	CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
516	IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[ 1,1,irec])
517	IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
518	CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
519	END DO
520
521	!--- Flip latitudes: ascending in input file, descending in "rlatu".
522	IF(n_dim==3) THEN
523	v(1,:,:,:) = v(1,nlat:1:-1,:,:)
524	v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1) !--- Duplication
525	ELSE
526	v(:,:,:,:) = v(:,nlat:1:-1,:,:)
527	END IF
528
529	END SUBROUTINE get_3Dfields
530	!
531	!-------------------------------------------------------------------------------
532
533
534
535	!-------------------------------------------------------------------------------
536	!
537	FUNCTION get_SigTrop(ih,it) RESULT(out)
538	!
539	!-------------------------------------------------------------------------------
540	! Arguments:
541	REAL :: out
542	INTEGER, INTENT(IN) :: ih
543	INTEGER, INTENT(OUT) :: it
544	!-------------------------------------------------------------------------------
545	!--- Pressure at tropopause read from the forcing file
546	IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
547
548	!--- Chemical tropopause definition based on a particular threshold
549	IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih))
550	ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))
551	ELSE ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF
552	out=out/Pgnd_in(ih)
553
554	END FUNCTION get_SigTrop
555	!
556	!-------------------------------------------------------------------------------
557
558
559	!-------------------------------------------------------------------------------
560	!
561	FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
562	!
563	!-------------------------------------------------------------------------------
564	! Purpose: Determine the tropopause using chemical criterium, ie the pressure
565	! at which the ozone concentration reaches a certain level.
566	!-------------------------------------------------------------------------------
567	! Remarks:
568	! * Input field is upside down (increasing pressure // increasing vertical idx)
569	! The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
570	! where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
571	! the O3 profile in this zone is decreasing with pressure.
572	! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
573	! determined using an empirical law roughly derived from ... & al.
574	!-------------------------------------------------------------------------------
575	! Arguments:
576	REAL :: out !--- Pressure at tropopause
577	INTEGER, INTENT(IN) :: ih !--- Horizontal index
578	INTEGER, INTENT(OUT) :: it !--- Index of tropopause layer
579	INTEGER, INTENT(IN) :: it0 !--- Idx: higher than tropopause
580	REAL, INTENT(IN) :: pres(:) !--- Pressure profile, increasing
581	REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV)
582	REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause
583	!-------------------------------------------------------------------------------
584	! Local variables:
585	REAL :: o3t !--- Ozone concent. at tropopause
586	REAL :: al !--- Interpolation coefficient
587	REAL :: coef !--- Coeff of latitude modulation
588	REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause
589	!-------------------------------------------------------------------------------
590	!--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
591	IF(PRESENT(o3trop)) THEN !=== THRESHOLD FROM ARGUMENTS
592	o3t=o3trop
593	ELSE !=== EMPIRICAL LAW
594	coef = TANH(lat_in(ih)/co3(3)) !--- Latitude modulation
595	coef = SIN (2.RPI(y_frac-1./6.)) * coef !--- Seasonal modulation
596	o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization
597	END IF
598
599	!--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
600	it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
601	al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
602	IF(Ploc=='C') out = pres(it)*(1.-al) pres(it+1)**al
603	IF(Ploc=='I') out = SQRT(pres(it)*(1.-al) pres(it+2)*al pres(it+1))
604	it = locate(pres(:), out) !--- pres(it)<Ptrop<pres(it+1)
605
606	END FUNCTION PTrop_chem
607	!
608	!-------------------------------------------------------------------------------
609
610
611	!-------------------------------------------------------------------------------
612	!
613	FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
614	!
615	!-------------------------------------------------------------------------------
616	IMPLICIT NONE
617	!-------------------------------------------------------------------------------
618	! Arguments:
619	LOGICAL :: out !--- .T. => some wrong values
620	REAL, INTENT(IN) :: o3col(:), lon, lat
621	INTEGER, INTENT(IN) :: ilev0
622	CHARACTER(LEN=*), INTENT(IN) :: layer
623	REAL, OPTIONAL, INTENT(IN) :: vmin
624	REAL, OPTIONAL, INTENT(IN) :: vmax
625	!-------------------------------------------------------------------------------
626	! Local variables:
627	INTEGER :: k
628	LOGICAL :: lmin, lmax
629	REAL :: fac
630	CHARACTER(LEN=6) :: unt
631	!-------------------------------------------------------------------------------
632	!--- NOTHING TO DO
633	lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
634	lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
635	out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
636
637	!--- SOME TOO LOW VALUES FOUND
638	IF(lmin) THEN
639	CALL unitc(vmin,unt,fac)
640	DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
641	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
642	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
643	fac*vmin,unt//' in '//TRIM(layer)
644	END DO
645	END IF
646
647	!--- SOME TOO HIGH VALUES FOUND
648	IF(lmax) THEN
649	CALL unitc(vmax,unt,fac)
650	DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
651	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
652	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
653	fac*vmax,unt//' in '//TRIM(layer)
654	END DO
655	END IF
656
657	END FUNCTION check_ozone
658	!
659	!-------------------------------------------------------------------------------
660
661
662	!-------------------------------------------------------------------------------
663	!
664	SUBROUTINE unitc(val,unt,fac)
665	!
666	!-------------------------------------------------------------------------------
667	IMPLICIT NONE
668	!-------------------------------------------------------------------------------
669	! Arguments:
670	REAL, INTENT(IN) :: val
671	CHARACTER(LEN=6), INTENT(OUT) :: unt
672	REAL, INTENT(OUT) :: fac
673	!-------------------------------------------------------------------------------
674	! Local variables:
675	INTEGER :: ndgt
676	!-------------------------------------------------------------------------------
677	ndgt=3*FLOOR(LOG10(val)/3.)
678	SELECT CASE(ndgt)
679	CASE(-6); unt=' ppmV '; fac=1.E6
680	CASE(-9); unt=' ppbV '; fac=1.E9
681	CASE DEFAULT; unt=' vmr '; fac=1.0
682	END SELECT
683
684	END SUBROUTINE unitc
685	!
686	!-------------------------------------------------------------------------------
687
688
689	!-------------------------------------------------------------------------------
690	!
691	SUBROUTINE msg(ll,str,sub)
692	!
693	!-------------------------------------------------------------------------------
694	USE print_control_mod, ONLY: lunout
695	IMPLICIT NONE
696	!-------------------------------------------------------------------------------
697	! Arguments:
698	LOGICAL, INTENT(IN) :: ll
699	CHARACTER(LEN=*), INTENT(IN) :: str
700	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
701	!-------------------------------------------------------------------------------
702	IF(.NOT.ll) RETURN
703	IF(PRESENT(sub)) THEN
704	WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
705	ELSE
706	WRITE(lunout,*)TRIM(str)
707	END IF
708
709	END SUBROUTINE msg
710	!
711	!-------------------------------------------------------------------------------
712
713	END SUBROUTINE regr_pr_time_av
714	!
715	!-------------------------------------------------------------------------------
716
717	END MODULE regr_pr_time_av_m

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

Download in other formats:

Original Format