Context Navigation

regr_pr_time_av_m.F90 @ 3435

Last change on this file since 3435 was 3435, checked in by Laurent Fairhead, 5 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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	END IF
252	!$OMP END MASTER
253	IF(lfirst) THEN
254	lfirst=.FALSE.; CALL bcast(lfirst)
255	IF(lAdjTro) CALL bcast(itrp0)
256	CALL bcast(lPrSfile); CALL bcast(lPrTfile)
257	CALL bcast(lO3Tfile); CALL bcast(linterp)
258	END IF
259
260	IF (is_master) THEN
261	ALLOCATE(ind_cell_glo_glo(klon_glo))
262	ELSE
263	ALLOCATE(ind_cell_glo_glo(0))
264	ENDIF
265	CALL gather(ind_cell_glo,ind_cell_glo_glo)
266	IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
267
268	CALL scatter2d(v1,v2)
269
270	!--- No "ps" in input file => assumed to be equal to current LMDZ ground press
271	IF(lPrSfile) THEN
272	IF (is_master .AND. grid_type==unstructured) pg1(:,:)=pg1(:,ind_cell_glo_glo(:))
273	CALL scatter2d(pg1,Pgnd_in)
274	ELSE
275	Pgnd_in=Pre_ou(:,1)
276	END IF
277
278	IF(lPrTfile) THEN
279	IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
280	CALL scatter2d(pt1,Ptrp_in)
281	ENDIF
282
283	IF(lO3Tfile) THEN
284	IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
285	CALL scatter2d(ot1,Otrp_in)
286	ENDIF
287	!--- No ground pressure in input file => choose it to be the one of LMDZ
288	IF(lAdjTro.AND..NOT.lPrSfile) Pgnd_in(:)=Pgrnd_ou(:)
289
290	!--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
291	IF(.NOT.lAdjTro) THEN
292	DO i=1,klon
293	Pres_ou=Pre_ou(i,SIZE(Pre_ou,2):1:-1) !--- pplay & paprs are decreasing
294	IF(Ploc=='C') CALL regr_lint (1,v2(i,:,:), LOG(Pre_in(:)), &
295	LOG(Pres_ou(:)), v3(i,nlev_ou:1:-1,:))
296	IF(Ploc=='I') CALL regr_conserv(1,v2(i,:,:), Pre_in(:) , &
297	Pres_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
298	END DO
299	!-------------------------------------------------------------------------------
300	ELSE !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
301	!-------------------------------------------------------------------------------
302	y_frac=(REAL(days_elapsed)+jH_cur)/year_len
303
304	!--- OUTPUT SIGMA LEVELS
305	DO i=1,klon
306
307	!--- INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN
308	Pres_ou = Pre_ou(i,SIZE(Pre_ou,2):1:-1)!--- pplay & paprs are decreasing
309	Sig_in(:) = Pre_in (:)/Pgnd_in(i) !--- increasing values
310	Sig_ou(:) = Pres_ou(:)/Pgnd_ou(i) !--- increasing values
311
312	!--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
313	! to keep tropopause pressure realistic ; high values are usually due to
314	! ozone hole fooling the crude chemical tropopause detection algorithm.
315	SigT_in = get_SigTrop(i,itrp)
316	SigT_in=MIN(SigT_in,ChemPTrMax/Pgnd_in(i)) !--- too low value filtered
317	SigT_in=MAX(SigT_in,ChemPTrMin/Pgnd_ou(i)) !--- too high value filtered
318
319	!--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
320	! dynamical tropopause (especially in filaments) are conterbalanced with
321	! a filter ensuring it stays within a certain distance around input (file)
322	! tropopause, hence avoiding avoid a too thick stretched region ; a final
323	! extra-safety filter keeps the tropopause pressure value realistic.
324	SigT_ou = Ptrp_ou(i)/Pgnd_ou(i)
325	IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low value w/r input
326	IF(SigT_ou>SigT_inrho) SigT_ou=SigT_inrho !--- too high value w/r input
327	SigT_ou=MIN(SigT_ou,DynPTrMax/Pgnd_ou(i)) !--- too low value filtered
328	SigT_ou=MAX(SigT_ou,DynPTrMin/Pgnd_ou(i)) !--- too high value filtered
329	Ptrop_ou(i)=SigT_ou*Pgnd_ou(i)
330	iout = locate(Sig_ou(:),SigT_ou)
331
332	!--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
333	alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
334
335	!--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS
336	Sig_bo0 = MAX(SigT_in,SigT_ou) !--- lowest tropopause
337	Sig_to0 = MIN(SigT_in,SigT_ou) !--- highest tropopause
338	beta = (Sig_bo0/Sig_to0)**gamm !--- stretching exponent
339	Sig_bot = MIN(Sig_bo0beta,0.1(9.+Sig_bo0)) !--- must be <1
340	ibot = locate(Sig_ou(:),Sig_bot) !--- layer index
341	IF(ibot-iout<2) THEN !--- at least one layer thick
342	ibot=MIN(iout+2,nlev_ou); Sig_bot=Sig_ou(ibot)
343	END IF
344	Sig_top = Sig_to0/beta !--- upper bound
345	itop = locate(Sig_ou(:),Sig_top) !--- layer index
346	IF(iout-itop<2) THEN !--- at least one layer thick
347	itop=MAX(iout-2,1); Sig_top=Sig_ou(itop)
348	END IF
349
350	!--- STRETCHING POWER LAW LOCALIZATION FUNCTION:
351	! 0 in [0,Sig_top] 0->1 in [Sig_top,SigT_ou]
352	! 0 in [Sig_bot,1] 1->0 in [SigT_ou, Sig_bot]
353	phi(:)=0.
354	phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&
355	*LOG(SigT_ou)/LOG(SigT_ou/Sig_top)
356	phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&
357	*LOG(SigT_ou)/LOG(SigT_ou/Sig_bot)
358
359	!--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
360	Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)*(alphaphi(:))
361
362	!--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
363	IF(Ploc=='C') CALL regr_lint (1, v2(i,:,:), LOG(Pre_in(:)), &
364	LOG(Pstr_ou(:)), v3(i,nlev_ou:1:-1,:))
365	IF(Ploc=='I') CALL regr_conserv(1, v2(i,:,:), Pre_in(:) , &
366	Pstr_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
367
368	!--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
369	i0=nlev_ou-locate(Pres_ou(:),Ptrop_ou(i))+1
370	ll=check_ozone(v3(i, 1:i0-1 ,1),lon_in(i),lat_in(i),1 ,'troposphere', &
371	5.E-9,3.0E-7)
372	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
373	ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
374	5.E-8,1.5E-5)
375	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
376
377	IF(ldebug) THEN
378	dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)
379	dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)
380	Ptrop_in (i) = SigT_in*Pgnd_in(i)
381	Pstrn_ou(i,:)= Pstr_ou
382	phii(i,:) = phi(:)
383	Ptrop_ef(i) = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp), &
384	Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)
385	END IF
386	END DO
387	END IF
388	IF(ldebug.AND.lAdjTro) THEN
389	CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres
390	CALL writefield_phy('dzStrain' ,dzStrain ,1) !--- Strained thickness
391	CALL writefield_phy('dzStrain0',dzStrain0,1) !--- Tropopauses distance
392	CALL writefield_phy('phi',phii,nlev_ou) !--- Localization function
393	!--- Tropopauses pressures:
394	CALL writefield_phy('PreTr_in',Ptrop_in,1) !--- Input and effective
395	CALL writefield_phy('PreTr_ou',Ptrop_ou,1) !--- LMDz dyn tropopause
396	CALL writefield_phy('PreTr_ef',Ptrop_ef,1) !--- Effective chem tropop
397	END IF
398	IF(ldebug) THEN
399	CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
400	CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field
401	CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure
402	END IF
403
404	CONTAINS
405
406
407	!-------------------------------------------------------------------------------
408	!
409	SUBROUTINE update_fields()
410	!
411	!-------------------------------------------------------------------------------
412	IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION
413	CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
414	irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number
415	!--- MIN -> Security in the unlikely case of roundup errors.
416	CALL get_3Dfields(v1) !--- Read ozone field(s)
417	IF(lAdjTro) THEN !--- Additional files for fields strain
418	IF(lPrSfile) CALL get_2Dfield(pg1,"ps")
419	IF(lPrTfile) CALL get_2Dfield(pt1,"tropopause_air_pressure")
420	IF(lO3Tfile) CALL get_2Dfield(ot1,"tro3_at_tropopause")
421	END IF
422	ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS
423	IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
424	IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
425	CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
426	IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS
427	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&
428	//'reading record ',irec,' (time ',time_in(irec),')'
429	CALL msg(.TRUE.,str,sub)
430	CALL get_3Dfields(v1m) !--- Read ozone field(s)
431	IF(lAdjTro) THEN !--- Additional files for fields strain
432	IF(lPrSfile) CALL get_2Dfield(pgm,"ps")
433	IF(lPrTfile) CALL get_2Dfield(ptm,"tropopause_air_pressure")
434	IF(lO3Tfile) CALL get_2Dfield(otm,"tro3_at_tropopause")
435	END IF
436	ELSE !=== SHIFT FIELDS
437	irec=irec+1
438	WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&
439	//' current next one (',irec,', time ',time_in(irec),')'
440	CALL msg(.TRUE.,str,sub)
441	v1m=v1p !--- Ozone fields
442	IF(lAdjTro) THEN !--- Additional files for fields strain
443	IF(lPrSfile) pgm=pgp !--- Surface pressure
444	IF(lPrTfile) ptm=ptp !--- Tropopause pressure
445	IF(lO3Tfile) otm=otp !--- Tropopause ozone
446	END IF
447	END IF
448	irec=irec+1
449	WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'&
450	,irec,' (time ',time_in(irec),')'
451	CALL msg(.TRUE.,str,sub)
452	CALL get_3Dfields(v1p) !--- Read ozone field(s)
453	IF(lAdjTro) THEN !--- Additional files for fields strain
454	IF(lPrSfile) CALL get_2Dfield(pgp,"ps")
455	IF(lPrTfile) CALL get_2Dfield(ptp,"tropopause_air_pressure")
456	IF(lO3Tfile) CALL get_2Dfield(otp,"tro3_at_tropopause")
457	END IF
458	irec=irec-1
459	END IF
460
461	END SUBROUTINE update_fields
462	!
463	!-------------------------------------------------------------------------------
464
465
466	!-------------------------------------------------------------------------------
467	!
468	SUBROUTINE get_2Dfield(v,var)
469	!
470	!-------------------------------------------------------------------------------
471	! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
472	! from the input file.
473	! Remark: In case the field is zonal, it is duplicated along first dimension.
474	!-------------------------------------------------------------------------------
475	! Arguments:
476	REAL, INTENT(INOUT) :: v(:,:)
477	CHARACTER(LEN=*), INTENT(IN) :: var
478	!-------------------------------------------------------------------------------
479	CALL NF95_INQ_VARID(fID, TRIM(var), vID)
480	CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
481	IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[ 1,irec])
482	IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
483	CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
484
485	!--- Flip latitudes: ascending in input file, descending in "rlatu".
486	IF(n_dim==3) THEN
487	v(1,:) = v(1,nlat:1:-1)
488	v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1) !--- Duplication
489	ELSE
490	v(:,:) = v(:,nlat:1:-1)
491	END IF
492
493	END SUBROUTINE get_2Dfield
494	!
495	!-------------------------------------------------------------------------------
496
497
498	!-------------------------------------------------------------------------------
499	!
500	SUBROUTINE get_3Dfields(v)
501	!
502	!-------------------------------------------------------------------------------
503	! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
504	! from the input file. Fields are stacked on fourth dimension.
505	! Remark: In case the field is zonal, it is duplicated along first dimension.
506	!-------------------------------------------------------------------------------
507	! Arguments:
508	REAL, INTENT(INOUT) :: v(:,:,:,:)
509	!-------------------------------------------------------------------------------
510	DO i=1,SIZE(nam)
511	CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
512	CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
513	IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[ 1,1,irec])
514	IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
515	CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
516	END DO
517
518	!--- Flip latitudes: ascending in input file, descending in "rlatu".
519	IF(n_dim==3) THEN
520	v(1,:,:,:) = v(1,nlat:1:-1,:,:)
521	v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1) !--- Duplication
522	ELSE
523	v(:,:,:,:) = v(:,nlat:1:-1,:,:)
524	END IF
525
526	END SUBROUTINE get_3Dfields
527	!
528	!-------------------------------------------------------------------------------
529
530
531
532	!-------------------------------------------------------------------------------
533	!
534	FUNCTION get_SigTrop(ih,it) RESULT(out)
535	!
536	!-------------------------------------------------------------------------------
537	! Arguments:
538	REAL :: out
539	INTEGER, INTENT(IN) :: ih
540	INTEGER, INTENT(OUT) :: it
541	!-------------------------------------------------------------------------------
542	!--- Pressure at tropopause read from the forcing file
543	IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
544
545	!--- Chemical tropopause definition based on a particular threshold
546	IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih))
547	ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))
548	ELSE ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF
549	out=out/Pgnd_in(ih)
550
551	END FUNCTION get_SigTrop
552	!
553	!-------------------------------------------------------------------------------
554
555
556	!-------------------------------------------------------------------------------
557	!
558	FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
559	!
560	!-------------------------------------------------------------------------------
561	! Purpose: Determine the tropopause using chemical criterium, ie the pressure
562	! at which the ozone concentration reaches a certain level.
563	!-------------------------------------------------------------------------------
564	! Remarks:
565	! * Input field is upside down (increasing pressure // increasing vertical idx)
566	! The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
567	! where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
568	! the O3 profile in this zone is decreasing with pressure.
569	! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
570	! determined using an empirical law roughly derived from ... & al.
571	!-------------------------------------------------------------------------------
572	! Arguments:
573	REAL :: out !--- Pressure at tropopause
574	INTEGER, INTENT(IN) :: ih !--- Horizontal index
575	INTEGER, INTENT(OUT) :: it !--- Index of tropopause layer
576	INTEGER, INTENT(IN) :: it0 !--- Idx: higher than tropopause
577	REAL, INTENT(IN) :: pres(:) !--- Pressure profile, increasing
578	REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV)
579	REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause
580	!-------------------------------------------------------------------------------
581	! Local variables:
582	REAL :: o3t !--- Ozone concent. at tropopause
583	REAL :: al !--- Interpolation coefficient
584	REAL :: coef !--- Coeff of latitude modulation
585	REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause
586	!-------------------------------------------------------------------------------
587	!--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
588	IF(PRESENT(o3trop)) THEN !=== THRESHOLD FROM ARGUMENTS
589	o3t=o3trop
590	ELSE !=== EMPIRICAL LAW
591	coef = TANH(lat_in(ih)/co3(3)) !--- Latitude modulation
592	coef = SIN (2.RPI(y_frac-1./6.)) * coef !--- Seasonal modulation
593	o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization
594	END IF
595
596	!--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
597	it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
598	al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
599	IF(Ploc=='C') out = pres(it)*(1.-al) pres(it+1)**al
600	IF(Ploc=='I') out = SQRT(pres(it)*(1.-al) pres(it+2)*al pres(it+1))
601	it = locate(pres(:), out) !--- pres(it)<Ptrop<pres(it+1)
602
603	END FUNCTION PTrop_chem
604	!
605	!-------------------------------------------------------------------------------
606
607
608	!-------------------------------------------------------------------------------
609	!
610	FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
611	!
612	!-------------------------------------------------------------------------------
613	IMPLICIT NONE
614	!-------------------------------------------------------------------------------
615	! Arguments:
616	LOGICAL :: out !--- .T. => some wrong values
617	REAL, INTENT(IN) :: o3col(:), lon, lat
618	INTEGER, INTENT(IN) :: ilev0
619	CHARACTER(LEN=*), INTENT(IN) :: layer
620	REAL, OPTIONAL, INTENT(IN) :: vmin
621	REAL, OPTIONAL, INTENT(IN) :: vmax
622	!-------------------------------------------------------------------------------
623	! Local variables:
624	INTEGER :: k
625	LOGICAL :: lmin, lmax
626	REAL :: fac
627	CHARACTER(LEN=6) :: unt
628	!-------------------------------------------------------------------------------
629	!--- NOTHING TO DO
630	lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
631	lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
632	out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
633
634	!--- SOME TOO LOW VALUES FOUND
635	IF(lmin) THEN
636	CALL unitc(vmin,unt,fac)
637	DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
638	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
639	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
640	fac*vmin,unt//' in '//TRIM(layer)
641	END DO
642	END IF
643
644	!--- SOME TOO HIGH VALUES FOUND
645	IF(lmax) THEN
646	CALL unitc(vmax,unt,fac)
647	DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
648	WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' &
649	//TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
650	fac*vmax,unt//' in '//TRIM(layer)
651	END DO
652	END IF
653
654	END FUNCTION check_ozone
655	!
656	!-------------------------------------------------------------------------------
657
658
659	!-------------------------------------------------------------------------------
660	!
661	SUBROUTINE unitc(val,unt,fac)
662	!
663	!-------------------------------------------------------------------------------
664	IMPLICIT NONE
665	!-------------------------------------------------------------------------------
666	! Arguments:
667	REAL, INTENT(IN) :: val
668	CHARACTER(LEN=6), INTENT(OUT) :: unt
669	REAL, INTENT(OUT) :: fac
670	!-------------------------------------------------------------------------------
671	! Local variables:
672	INTEGER :: ndgt
673	!-------------------------------------------------------------------------------
674	ndgt=3*FLOOR(LOG10(val)/3.)
675	SELECT CASE(ndgt)
676	CASE(-6); unt=' ppmV '; fac=1.E6
677	CASE(-9); unt=' ppbV '; fac=1.E9
678	CASE DEFAULT; unt=' vmr '; fac=1.0
679	END SELECT
680
681	END SUBROUTINE unitc
682	!
683	!-------------------------------------------------------------------------------
684
685
686	!-------------------------------------------------------------------------------
687	!
688	SUBROUTINE msg(ll,str,sub)
689	!
690	!-------------------------------------------------------------------------------
691	USE print_control_mod, ONLY: lunout
692	IMPLICIT NONE
693	!-------------------------------------------------------------------------------
694	! Arguments:
695	LOGICAL, INTENT(IN) :: ll
696	CHARACTER(LEN=*), INTENT(IN) :: str
697	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
698	!-------------------------------------------------------------------------------
699	IF(.NOT.ll) RETURN
700	IF(PRESENT(sub)) THEN
701	WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
702	ELSE
703	WRITE(lunout,*)TRIM(str)
704	END IF
705
706	END SUBROUTINE msg
707	!
708	!-------------------------------------------------------------------------------
709
710	END SUBROUTINE regr_pr_time_av
711	!
712	!-------------------------------------------------------------------------------
713
714	END MODULE regr_pr_time_av_m

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

Download in other formats:

Original Format