Context Navigation

regr_pr_time_av_m.f90 @ 5319

Last change on this file since 5319 was 5285, checked in by abarral, 10 days ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 34.6 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, &
116	NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, nf95_get_var
117	USE netcdf, ONLY: NF90_INQ_VARID, 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	USE yomcst_mod_h
131	USE clesphys_mod_h
132	!-------------------------------------------------------------------------------
133	! Arguments:
134	INTEGER, INTENT(IN) :: fID !--- NetCDF file ID
135	CHARACTER(LEN=13), INTENT(IN) :: nam(:) !--- NetCDF variables names
136	REAL, INTENT(IN) :: julien !--- Days since Jan 1st
137	CHARACTER(LEN=1), INTENT(IN) :: Ploc !--- Pressures locations
138	!--- File/LMDZ (resp. decreasing & increasing order) pressure, Pa
139	! At cells centers or interfaces depending on "Ploc" keyword (C/I)
140	REAL, INTENT(IN) :: Pre_in(:) !--- in: file (nlev_in[+1])
141	REAL, INTENT(IN) :: Pre_ou(:,:) !--- out: LMDZ (klon,nlev_ou[+1])
142	REAL, INTENT(OUT) :: v3(:,:,:) !--- Regr. fld (klon,nlev_ou,n_var)
143	REAL, INTENT(IN), OPTIONAL :: Pgnd_ou(:) !--- LMDZ ground pressure (klon)
144	REAL, INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
145	! since Jan 1 of current year
146	REAL, INTENT(IN), OPTIONAL :: lon_in(:) !--- File longitudes vector (klon)
147	REAL, INTENT(IN), OPTIONAL :: lat_in(:) !--- File latitudes vector (klon)
148	REAL, INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres (klon)
149	!-------------------------------------------------------------------------------
150	! Local variables:
151
152	CHARACTER(LEN=80) :: sub
153	CHARACTER(LEN=320) :: str
154	INTEGER :: vID, ncerr, n_var, ibot, iout, nn
155	INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
156	LOGICAL :: lAdjTro !--- Need to adjust tropopause
157	REAL :: y_frac !--- Elapsed year fraction
158	REAL :: alpha, beta, al !--- For stretching/interpolation
159	REAL :: SigT_in, SigT_ou !--- Input and output tropopauses
160	REAL :: Sig_bot, Sig_top !--- Bounds of quasi-hat function
161	REAL :: Sig_bo0, Sig_to0 !--- Lower/upper tropopauses
162	REAL :: Sig_in (SIZE(Pre_in)) !--- Input field sigma levels
163	REAL :: Sig_ou (SIZE(Pre_ou,2)) !--- Output LMDZ sigma levels
164	REAL :: phi (SIZE(Pre_ou,2)) !--- Stretching exponent anomaly
165	REAL :: Pstr_ou(SIZE(Pre_ou,2)) !--- Stretched pressure levels
166	REAL :: Pres_ou(SIZE(Pre_ou,2)) !--- Pre_ou(i,:), reversed order
167	REAL, DIMENSION(nlon, nlat) :: pg1, pt1, ot1
168	REAL, DIMENSION(klon) :: Pgnd_in, Ptrp_in, Otrp_in
169	REAL, DIMENSION(klon) :: Ptrop_ou, Pgrnd_ou
170	! * The following fields are scattered to the partial "physics" horizontal grid.
171	REAL, POINTER :: v2(:,:,:) !--- Current time ozone fields
172	! In the 2D file case, the values are the same for all longitudes.
173	! "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)" and name "nam(l)"
174	! Both are: * if Ploc=='I' in pressure interval "press_in_edg(k:k+1)"
175	! * if Ploc=='C' at pressure "press_in_cen(k)"
176	REAL, TARGET :: &
177	v2i(klon,SIZE(Pre_in)-1,SIZE(nam)), & !--- v2 in Ploc=='I' case
178	v2c(klon,SIZE(Pre_in) ,SIZE(nam)) !--- v2 in Ploc=='C' case
179	INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
180	LOGICAL :: ll
181	!--- For debug
182	REAL, DIMENSION(klon) :: Ptrop_in, Ptrop_ef
183	REAL, DIMENSION(klon) :: dzStrain, dzStrain0
184	REAL, DIMENSION(klon,SIZE(Pre_ou,2)) :: Pstrn_ou, phii
185	!-------------------------------------------------------------------------------
186	sub="regr_pr_time_av"
187	nlev_in=SIZE(Pre_in); IF(Ploc=='I') nlev_in=nlev_in-1
188	IF(Ploc=='I') THEN; v2 => v2i; ELSE; v2 => v2c; END IF
189	CALL assert(SIZE(Pre_ou,1)==klon,TRIM(sub)//" Pre_ou klon")
190	CALL assert(SIZE(v3,1)==klon, TRIM(sub)//" v3 klon")
191	CALL assert(SIZE(v3,2)==nlev_ou, TRIM(sub)//" v3 nlev_ou")
192	IF(Ploc=='I') CALL assert(SIZE(Pre_ou,2)==nlev_ou+1,TRIM(sub)//" Pre_ou nlev_ou+1")
193	IF(Ploc=='C') CALL assert(SIZE(Pre_ou,2)==nlev_ou ,TRIM(sub)//" Pre_ou nlev_ou")
194	n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
195	IF(PRESENT(Pgnd_ou)) CALL assert(SIZE(Pgnd_ou)==klon,TRIM(sub)//" Pgnd_ou klon")
196	IF(PRESENT(lon_in)) CALL assert(SIZE(lon_in )==klon,TRIM(sub)//" lon_in klon")
197	IF(PRESENT(lat_in)) CALL assert(SIZE(lat_in )==klon,TRIM(sub)//" lat_in klon")
198	IF(PRESENT(Ptrp_ou)) CALL assert(SIZE(Ptrp_ou)==klon,TRIM(sub)//" Ptrp_ou klon")
199	lAdjTro=PRESENT(Ptrp_ou)
200	IF(lAdjTro) THEN
201	IF(.NOT.PRESENT(lat_in)) &
202	CALL abort_physic(sub, 'Missing lat_in (required if adjust_tropopause=T)', 1)
203	IF(.NOT.PRESENT(Pgnd_ou).AND.Ploc=='C') &
204	CALL abort_physic(sub, 'Missing ground Pr(required if adjust_tropopause=T)', 1)
205	IF(PRESENT(Pgnd_ou)) THEN; Pgrnd_ou=Pgnd_ou; ELSE; Pgrnd_ou=Pre_ou(:,1); END IF
206	END IF
207
208	!$OMP MASTER
209	IF(is_mpi_root) THEN
210
211	!=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
212	IF(lfirst) THEN
213	lPrSfile=lAdjTro.AND.NF90_INQ_VARID(fID,"ps" ,vID)==NF90_NOERR
214	lPrTfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
215	lO3Tfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tro3_at_tropopause" ,vID)==NF90_NOERR
216	CALL NF95_INQ_DIMID(fID,"time",vID)
217	CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
218	linterp=PRESENT(time_in).AND.ntim_in==14
219	ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
220	IF(linterp) THEN
221	ALLOCATE(v1m(nlon,nlat,nlev_in,n_var),v1p(nlon,nlat,nlev_in,n_var))
222	IF(lPrSfile) ALLOCATE(pgm(nlon,nlat),pgp(nlon,nlat))
223	IF(lPrTfile) ALLOCATE(ptm(nlon,nlat),ptp(nlon,nlat))
224	IF(lO3Tfile) ALLOCATE(otm(nlon,nlat),otp(nlon,nlat))
225	END IF
226	!--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
227	IF(lAdjTro) itrp0=locate(Pre_in,pTropUp)
228	CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.' ,sub)
229	CALL msg(lPrSfile,'Using GROUND PRESSURE from input O3 forcing file.',sub)
230	CALL msg(lAdjTro ,'o3 forcing file tropopause location uses:' ,sub)
231	IF(lPrTfile) THEN; str=' INPUT FILE PRESSURE'
232	ELSE IF(lO3Tfile) THEN; str=' INPUT FILE O3 CONCENTRATION'
233	ELSE IF(lO3Tpara) THEN; str=' PARAMETRIZED O3 concentration'
234	ELSE; str=' CONSTANT O3 concentration'; END IF
235	CALL msg(lAdjTro,TRIM(str)//' at tropopause')
236	END IF
237
238	!=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
239	CALL update_fields()
240
241	!=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
242	IF(linterp) THEN
243	WRITE(str,'(a,f12.8,2(a,f5.1))')'Interpolating O3 at julian day ',julien,&
244	' from fields at times ',time_in(irec),' and ', time_in(irec+1)
245	CALL msg(.TRUE.,str,sub)
246	al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
247	v1=alv1m+(1.-al)v1p
248	IF(lPrSfile) pg1=alpgm+(1.-al)pgp
249	IF(lPrTfile) pt1=alptm+(1.-al)ptp
250	IF(lO3Tfile) ot1=alotm+(1.-al)otp
251	END IF
252	ELSE
253	IF (lfirst) ALLOCATE(v1(0,0,0,0))
254	END IF
255	!$OMP END MASTER
256	!$OMP BARRIER
257	IF(lfirst) THEN
258	lfirst=.FALSE.; CALL bcast(lfirst)
259	IF(lAdjTro) CALL bcast(itrp0)
260	CALL bcast(lPrSfile); CALL bcast(lPrTfile)
261	CALL bcast(lO3Tfile); CALL bcast(linterp)
262	END IF
263
264	IF (is_master) THEN
265	ALLOCATE(ind_cell_glo_glo(klon_glo))
266	ELSE
267	ALLOCATE(ind_cell_glo_glo(0))
268	ENDIF
269	CALL gather(ind_cell_glo,ind_cell_glo_glo)
270	IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
271
272	CALL scatter2d(v1,v2)
273
274	!--- No "ps" in input file => assumed to be equal to current LMDZ ground press
275	IF(lPrSfile) THEN
276	IF (is_master .AND. grid_type==unstructured) pg1(:,:)=pg1(:,ind_cell_glo_glo(:))
277	CALL scatter2d(pg1,Pgnd_in)
278	ELSE
279	Pgnd_in=Pre_ou(:,1)
280	END IF
281
282	IF(lPrTfile) THEN
283	IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
284	CALL scatter2d(pt1,Ptrp_in)
285	ENDIF
286
287	IF(lO3Tfile) THEN
288	IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
289	CALL scatter2d(ot1,Otrp_in)
290	ENDIF
291	!--- No ground pressure in input file => choose it to be the one of LMDZ
292	IF(lAdjTro.AND..NOT.lPrSfile) Pgnd_in(:)=Pgrnd_ou(:)
293
294	!--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
295	IF(.NOT.lAdjTro) THEN
296	DO i=1,klon
297	Pres_ou=Pre_ou(i,SIZE(Pre_ou,2):1:-1) !--- pplay & paprs are decreasing
298	IF(Ploc=='C') CALL regr_lint (1,v2(i,:,:), LOG(Pre_in(:)), &
299	LOG(Pres_ou(:)), v3(i,nlev_ou:1:-1,:))
300	IF(Ploc=='I') CALL regr_conserv(1,v2(i,:,:), Pre_in(:) , &
301	Pres_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
302	END DO
303	!-------------------------------------------------------------------------------
304	ELSE !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
305	!-------------------------------------------------------------------------------
306	y_frac=(REAL(days_elapsed)+jH_cur)/year_len
307
308	!--- OUTPUT SIGMA LEVELS
309	DO i=1,klon
310
311	!--- INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN
312	Pres_ou = Pre_ou(i,SIZE(Pre_ou,2):1:-1)!--- pplay & paprs are decreasing
313	Sig_in(:) = Pre_in (:)/Pgnd_in(i) !--- increasing values
314	Sig_ou(:) = Pres_ou(:)/Pgnd_ou(i) !--- increasing values
315
316	!--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
317	! to keep tropopause pressure realistic ; high values are usually due to
318	! ozone hole fooling the crude chemical tropopause detection algorithm.
319	SigT_in = get_SigTrop(i,itrp)
320	SigT_in=MIN(SigT_in,ChemPTrMax/Pgnd_in(i)) !--- too low value filtered
321	SigT_in=MAX(SigT_in,ChemPTrMin/Pgnd_ou(i)) !--- too high value filtered
322
323	!--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
324	! dynamical tropopause (especially in filaments) are conterbalanced with
325	! a filter ensuring it stays within a certain distance around input (file)
326	! tropopause, hence avoiding avoid a too thick stretched region ; a final
327	! extra-safety filter keeps the tropopause pressure value realistic.
328	SigT_ou = Ptrp_ou(i)/Pgnd_ou(i)
329	IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low value w/r input
330	IF(SigT_ou>SigT_inrho) SigT_ou=SigT_inrho !--- too high value w/r input
331	SigT_ou=MIN(SigT_ou,DynPTrMax/Pgnd_ou(i)) !--- too low value filtered
332	SigT_ou=MAX(SigT_ou,DynPTrMin/Pgnd_ou(i)) !--- too high value filtered
333	Ptrop_ou(i)=SigT_ou*Pgnd_ou(i)
334	iout = locate(Sig_ou(:),SigT_ou)
335
336	!--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
337	alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
338
339	!--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS
340	Sig_bo0 = MAX(SigT_in,SigT_ou) !--- lowest tropopause
341	Sig_to0 = MIN(SigT_in,SigT_ou) !--- highest tropopause
342	beta = (Sig_bo0/Sig_to0)**gamm !--- stretching exponent
343	Sig_bot = MIN(Sig_bo0beta,0.1(9.+Sig_bo0)) !--- must be <1
344	ibot = locate(Sig_ou(:),Sig_bot) !--- layer index
345	IF(ibot-iout<2) THEN !--- at least one layer thick
346	ibot=MIN(iout+2,nlev_ou); Sig_bot=Sig_ou(ibot)
347	END IF
348	Sig_top = Sig_to0/beta !--- upper bound
349	itop = locate(Sig_ou(:),Sig_top) !--- layer index
350	IF(iout-itop<2) THEN !--- at least one layer thick
351	itop=MAX(iout-2,1); Sig_top=Sig_ou(itop)
352	END IF
353
354	!--- STRETCHING POWER LAW LOCALIZATION FUNCTION:
355	! 0 in [0,Sig_top] 0->1 in [Sig_top,SigT_ou]
356	! 0 in [Sig_bot,1] 1->0 in [SigT_ou, Sig_bot]
357	phi(:)=0.
358	phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&
359	*LOG(SigT_ou)/LOG(SigT_ou/Sig_top)
360	phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&
361	*LOG(SigT_ou)/LOG(SigT_ou/Sig_bot)
362
363	!--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
364	Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)*(alphaphi(:))
365
366	!--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
367	IF(Ploc=='C') CALL regr_lint (1, v2(i,:,:), LOG(Pre_in(:)), &
368	LOG(Pstr_ou(:)), v3(i,nlev_ou:1:-1,:))
369	IF(Ploc=='I') CALL regr_conserv(1, v2(i,:,:), Pre_in(:) , &
370	Pstr_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
371
372	!--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
373	i0=nlev_ou-locate(Pres_ou(:),Ptrop_ou(i))+1
374	ll=check_ozone(v3(i, 1:i0-1 ,1),lon_in(i),lat_in(i),1 ,'troposphere', &
375	5.E-9,3.0E-7)
376	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
377	ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
378	5.E-8,1.5E-5)
379	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
380
381	IF(ldebug) THEN
382	dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)
383	dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)
384	Ptrop_in (i) = SigT_in*Pgnd_in(i)
385	Pstrn_ou(i,:)= Pstr_ou
386	phii(i,:) = phi(:)
387	Ptrop_ef(i) = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp), &
388	Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)
389	END IF
390	END DO
391	END IF
392	IF(ldebug.AND.lAdjTro) THEN
393	CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres
394	CALL writefield_phy('dzStrain' ,dzStrain ,1) !--- Strained thickness
395	CALL writefield_phy('dzStrain0',dzStrain0,1) !--- Tropopauses distance
396	CALL writefield_phy('phi',phii,nlev_ou) !--- Localization function
397	!--- Tropopauses pressures:
398	CALL writefield_phy('PreTr_in',Ptrop_in,1) !--- Input and effective
399	CALL writefield_phy('PreTr_ou',Ptrop_ou,1) !--- LMDz dyn tropopause
400	CALL writefield_phy('PreTr_ef',Ptrop_ef,1) !--- Effective chem tropop
401	END IF
402	IF(ldebug) THEN
403	CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
404	CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field
405	CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure
406	END IF
407
408	CONTAINS
409
410
411	!-------------------------------------------------------------------------------
412	!
413	SUBROUTINE update_fields()
414	!
415	!-------------------------------------------------------------------------------
416	IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION
417	CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
418	irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number
419	!--- MIN -> Security in the unlikely case of roundup errors.
420	CALL get_3Dfields(v1) !--- Read ozone field(s)
421	IF(lAdjTro) THEN !--- Additional files for fields strain
422	IF(lPrSfile) CALL get_2Dfield(pg1,"ps")
423	IF(lPrTfile) CALL get_2Dfield(pt1,"tropopause_air_pressure")
424	IF(lO3Tfile) CALL get_2Dfield(ot1,"tro3_at_tropopause")
425	END IF
426	ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS
427	IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
428	IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
429	CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
430	IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS
431	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&
432	//'reading record ',irec,' (time ',time_in(irec),')'
433	CALL msg(.TRUE.,str,sub)
434	CALL get_3Dfields(v1m) !--- Read ozone field(s)
435	IF(lAdjTro) THEN !--- Additional files for fields strain
436	IF(lPrSfile) CALL get_2Dfield(pgm,"ps")
437	IF(lPrTfile) CALL get_2Dfield(ptm,"tropopause_air_pressure")
438	IF(lO3Tfile) CALL get_2Dfield(otm,"tro3_at_tropopause")
439	END IF
440	ELSE !=== SHIFT FIELDS
441	irec=irec+1
442	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&
443	//' current next one (',irec,', time ',time_in(irec),')'
444	CALL msg(.TRUE.,str,sub)
445	v1m=v1p !--- Ozone fields
446	IF(lAdjTro) THEN !--- Additional files for fields strain
447	IF(lPrSfile) pgm=pgp !--- Surface pressure
448	IF(lPrTfile) ptm=ptp !--- Tropopause pressure
449	IF(lO3Tfile) otm=otp !--- Tropopause ozone
450	END IF
451	END IF
452	irec=irec+1
453	WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'&
454	,irec,' (time ',time_in(irec),')'
455	CALL msg(.TRUE.,str,sub)
456	CALL get_3Dfields(v1p) !--- Read ozone field(s)
457	IF(lAdjTro) THEN !--- Additional files for fields strain
458	IF(lPrSfile) CALL get_2Dfield(pgp,"ps")
459	IF(lPrTfile) CALL get_2Dfield(ptp,"tropopause_air_pressure")
460	IF(lO3Tfile) CALL get_2Dfield(otp,"tro3_at_tropopause")
461	END IF
462	irec=irec-1
463	END IF
464
465	END SUBROUTINE update_fields
466	!
467	!-------------------------------------------------------------------------------
468
469
470	!-------------------------------------------------------------------------------
471	!
472	SUBROUTINE get_2Dfield(v,var)
473	!
474	!-------------------------------------------------------------------------------
475	! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
476	! from the input file.
477	! Remark: In case the field is zonal, it is duplicated along first dimension.
478	!-------------------------------------------------------------------------------
479	! Arguments:
480	REAL, INTENT(INOUT) :: v(:,:)
481	CHARACTER(LEN=*), INTENT(IN) :: var
482	!-------------------------------------------------------------------------------
483	CALL NF95_INQ_VARID(fID, TRIM(var), vID)
484	CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
485	IF(n_dim==2) call NF95_GET_VAR(fID,vID,v(1,:), start=[ 1,irec])
486	IF(n_dim==3) call NF95_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
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) call NF95_GET_VAR(fID,vID,v(1,:,:,i), start=[ 1,1,irec])
517	IF(n_dim==4) call NF95_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
518	END DO
519
520	!--- Flip latitudes: ascending in input file, descending in "rlatu".
521	IF(n_dim==3) THEN
522	v(1,:,:,:) = v(1,nlat:1:-1,:,:)
523	v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1) !--- Duplication
524	ELSE
525	v(:,:,:,:) = v(:,nlat:1:-1,:,:)
526	END IF
527
528	END SUBROUTINE get_3Dfields
529	!
530	!-------------------------------------------------------------------------------
531
532
533
534	!-------------------------------------------------------------------------------
535	!
536	FUNCTION get_SigTrop(ih,it) RESULT(out)
537	!
538	!-------------------------------------------------------------------------------
539	! Arguments:
540	REAL :: out
541	INTEGER, INTENT(IN) :: ih
542	INTEGER, INTENT(OUT) :: it
543	!-------------------------------------------------------------------------------
544	!--- Pressure at tropopause read from the forcing file
545	IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
546
547	!--- Chemical tropopause definition based on a particular threshold
548	IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih))
549	ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))
550	ELSE ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF
551	out=out/Pgnd_in(ih)
552
553	END FUNCTION get_SigTrop
554	!
555	!-------------------------------------------------------------------------------
556
557
558	!-------------------------------------------------------------------------------
559	!
560	FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
561	!
562	!-------------------------------------------------------------------------------
563	! Purpose: Determine the tropopause using chemical criterium, ie the pressure
564	! at which the ozone concentration reaches a certain level.
565	!-------------------------------------------------------------------------------
566	! Remarks:
567	! * Input field is upside down (increasing pressure // increasing vertical idx)
568	! The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
569	! where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
570	! the O3 profile in this zone is decreasing with pressure.
571	! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
572	! determined using an empirical law roughly derived from ... & al.
573	!-------------------------------------------------------------------------------
574	! Arguments:
575	REAL :: out !--- Pressure at tropopause
576	INTEGER, INTENT(IN) :: ih !--- Horizontal index
577	INTEGER, INTENT(OUT) :: it !--- Index of tropopause layer
578	INTEGER, INTENT(IN) :: it0 !--- Idx: higher than tropopause
579	REAL, INTENT(IN) :: pres(:) !--- Pressure profile, increasing
580	REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV)
581	REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause
582	!-------------------------------------------------------------------------------
583	! Local variables:
584	REAL :: o3t !--- Ozone concent. at tropopause
585	REAL :: al !--- Interpolation coefficient
586	REAL :: coef !--- Coeff of latitude modulation
587	REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause
588	!-------------------------------------------------------------------------------
589	!--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
590	IF(PRESENT(o3trop)) THEN !=== THRESHOLD FROM ARGUMENTS
591	o3t=o3trop
592	ELSE !=== EMPIRICAL LAW
593	coef = TANH(lat_in(ih)/co3(3)) !--- Latitude modulation
594	coef = SIN (2.RPI(y_frac-1./6.)) * coef !--- Seasonal modulation
595	o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization
596	END IF
597
598	!--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
599	it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
600	al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
601	IF(Ploc=='C') out = pres(it)*(1.-al) pres(it+1)**al
602	IF(Ploc=='I') out = SQRT(pres(it)*(1.-al) pres(it+2)*al pres(it+1))
603	it = locate(pres(:), out) !--- pres(it)<Ptrop<pres(it+1)
604
605	END FUNCTION PTrop_chem
606	!
607	!-------------------------------------------------------------------------------
608
609
610	!-------------------------------------------------------------------------------
611	!
612	FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
613	!
614	!-------------------------------------------------------------------------------
615	IMPLICIT NONE
616	!-------------------------------------------------------------------------------
617	! Arguments:
618	LOGICAL :: out !--- .T. => some wrong values
619	REAL, INTENT(IN) :: o3col(:), lon, lat
620	INTEGER, INTENT(IN) :: ilev0
621	CHARACTER(LEN=*), INTENT(IN) :: layer
622	REAL, OPTIONAL, INTENT(IN) :: vmin
623	REAL, OPTIONAL, INTENT(IN) :: vmax
624	!-------------------------------------------------------------------------------
625	! Local variables:
626	INTEGER :: k
627	LOGICAL :: lmin, lmax
628	REAL :: fac
629	CHARACTER(LEN=6) :: unt
630	!-------------------------------------------------------------------------------
631	!--- NOTHING TO DO
632	lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
633	lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
634	out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
635
636	!--- SOME TOO LOW VALUES FOUND
637	IF(lmin) THEN
638	CALL unitc(vmin,unt,fac)
639	DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
640	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
641	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
642	fac*vmin,unt//' in '//TRIM(layer)
643	END DO
644	END IF
645
646	!--- SOME TOO HIGH VALUES FOUND
647	IF(lmax) THEN
648	CALL unitc(vmax,unt,fac)
649	DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
650	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
651	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
652	fac*vmax,unt//' in '//TRIM(layer)
653	END DO
654	END IF
655
656	END FUNCTION check_ozone
657	!
658	!-------------------------------------------------------------------------------
659
660
661	!-------------------------------------------------------------------------------
662	!
663	SUBROUTINE unitc(val,unt,fac)
664	!
665	!-------------------------------------------------------------------------------
666	IMPLICIT NONE
667	!-------------------------------------------------------------------------------
668	! Arguments:
669	REAL, INTENT(IN) :: val
670	CHARACTER(LEN=6), INTENT(OUT) :: unt
671	REAL, INTENT(OUT) :: fac
672	!-------------------------------------------------------------------------------
673	! Local variables:
674	INTEGER :: ndgt
675	!-------------------------------------------------------------------------------
676	ndgt=3*FLOOR(LOG10(val)/3.)
677	SELECT CASE(ndgt)
678	CASE(-6); unt=' ppmV '; fac=1.E6
679	CASE(-9); unt=' ppbV '; fac=1.E9
680	CASE DEFAULT; unt=' vmr '; fac=1.0
681	END SELECT
682
683	END SUBROUTINE unitc
684	!
685	!-------------------------------------------------------------------------------
686
687
688	!-------------------------------------------------------------------------------
689	!
690	SUBROUTINE msg(ll,str,sub)
691	!
692	!-------------------------------------------------------------------------------
693	USE print_control_mod, ONLY: lunout
694	IMPLICIT NONE
695	!-------------------------------------------------------------------------------
696	! Arguments:
697	LOGICAL, INTENT(IN) :: ll
698	CHARACTER(LEN=*), INTENT(IN) :: str
699	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
700	!-------------------------------------------------------------------------------
701	IF(.NOT.ll) RETURN
702	IF(PRESENT(sub)) THEN
703	WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
704	ELSE
705	WRITE(lunout,*)TRIM(str)
706	END IF
707
708	END SUBROUTINE msg
709	!
710	!-------------------------------------------------------------------------------
711
712	END SUBROUTINE regr_pr_time_av
713	!
714	!-------------------------------------------------------------------------------
715
716	END MODULE regr_pr_time_av_m

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

Download in other formats:

Original Format