Context Navigation

regr_pr_time_av_m.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
File size: 34.7 KB

Rev	Line
[2788]	1	! $Id$
	2	MODULE regr_pr_time_av_m
	3
	4	USE write_field_phy
	5	USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2968]	6	USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
[3069]	7	USE print_control_mod, ONLY: prt_level
[2788]	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.
[2968]	17	! * the forcing fields are stretched if the following arguments are present:
[3086]	18	! - "lat_in": input file latitudes.
	19	! - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
[2968]	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.
[2788]	22	!-------------------------------------------------------------------------------
	23	! Initial routine: regr_pr_av_m module (L. Guez).
	24	! Present version: David Cugnet ; corresponding additions:
	25	! - time interpolation
	26	! - 3D compliant
[2968]	27	! - vertical stretching of the field to allow tropopause matching between
	28	! input field and current lmdz field.
[2788]	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.
[3086]	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.
[2788]	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.
[2968]	48	! 2) computed using "tro3" and "tro3_at_tropopause' (if available).
[3086]	49	! 3) computed using "tro3" and a fixed threshold otherwise, constant or
	50	! determined using an empirical three parameters law:
[2968]	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.
[2788]	54	! * If available, the field "ps" (input file pressure at surface) is used.
	55	! If not, the current LMDZ ground pressure is taken instead.
[2968]	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.
[3086]	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)]".
[2968]	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
[3086]	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)
[2968]	66	! The stretching function has the following form:
[3086]	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.
[2968]	75	!
[2788]	76	! * The following fields are on the global "dynamics" grid, as read from files:
[3086]	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
[2788]	82	REAL, SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
	83	REAL, SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
[2968]	84	INTEGER, SAVE :: ntim_in !--- Records nb in input file
	85	INTEGER, SAVE :: itrp0 !--- idx above chem tropop.
[2788]	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
[3086]	90	LOGICAL, SAVE :: lPrSfile !--- Surface pressure flag
	91	LOGICAL, SAVE :: lPrTfile !--- Tropopause pressure flag
	92	LOGICAL, SAVE :: lO3Tfile !--- Tropopause ozone flag
[2788]	93	LOGICAL, SAVE :: lfirst=.TRUE. !--- First call flag
[2819]	94	!$OMP THREADPRIVATE(lfirst)
[3069]	95	REAL, PARAMETER :: pTropUp=9.E+3 !--- Value < tropopause pressure (Pa)
[3086]	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
[3087]	100	LOGICAL, PARAMETER :: ldebug=.FALSE.!--- Force writefield_phy multiple outputs
[3069]	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).
[2788]	105
	106	CONTAINS
	107
	108	!-------------------------------------------------------------------------------
	109	!
[3086]	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)
[2788]	112	!
	113	!-------------------------------------------------------------------------------
	114	USE dimphy, ONLY: klon
[2968]	115	USE netcdf95, ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
	116	NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
[2788]	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
[3605]	120	!! USE comvert_mod, ONLY: scaleheight
[2788]	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
[3605]	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
[2788]	128	USE phys_cal_mod, ONLY: calend, year_len, days_elapsed, jH_cur
[3605]	129	USE geometry_mod, ONLY: ind_cell_glo
[2788]	130	!-------------------------------------------------------------------------------
	131	! Arguments:
[3086]	132	INTEGER, INTENT(IN) :: fID !--- NetCDF file ID
[2788]	133	CHARACTER(LEN=13), INTENT(IN) :: nam(:) !--- NetCDF variables names
[3086]	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)
[2788]	142	REAL, INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
	143	! since Jan 1 of current year
[3069]	144	REAL, INTENT(IN), OPTIONAL :: lon_in(:) !--- File longitudes vector (klon)
	145	REAL, INTENT(IN), OPTIONAL :: lat_in(:) !--- File latitudes vector (klon)
[3086]	146	REAL, INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres (klon)
[2788]	147	!-------------------------------------------------------------------------------
	148	! Local variables:
	149	include "clesphys.h"
	150	include "YOMCST.h"
	151	CHARACTER(LEN=80) :: sub
[2968]	152	CHARACTER(LEN=320) :: str
[3086]	153	INTEGER :: vID, ncerr, n_var, ibot, iout, nn
	154	INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
[2820]	155	LOGICAL :: lAdjTro !--- Need to adjust tropopause
[2788]	156	REAL :: y_frac !--- Elapsed year fraction
[2968]	157	REAL :: alpha, beta, al !--- For stretching/interpolation
	158	REAL :: SigT_in, SigT_ou !--- Input and output tropopauses
[3086]	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
[2788]	171	! In the 2D file case, the values are the same for all longitudes.
[3086]	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
[3605]	178	INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
[2968]	179	LOGICAL :: ll
[3069]	180	!--- For debug
[3086]	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
[2788]	184	!-------------------------------------------------------------------------------
	185	sub="regr_pr_time_av"
[3086]	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")
[2968]	193	n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
[3086]	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
[2788]	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
[3086]	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
[2968]	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
[3086]	218	ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
[2788]	219	IF(linterp) THEN
[3086]	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))
[2788]	224	END IF
[2819]	225	!--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
[3086]	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')
[2788]	235	END IF
	236
	237	!=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
[2820]	238	CALL update_fields()
[2788]	239
	240	!=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
	241	IF(linterp) THEN
[3086]	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)
[2968]	244	CALL msg(.TRUE.,str,sub)
[2788]	245	al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
	246	v1=alv1m+(1.-al)v1p
[3086]	247	IF(lPrSfile) pg1=alpgm+(1.-al)pgp
	248	IF(lPrTfile) pt1=alptm+(1.-al)ptp
	249	IF(lO3Tfile) ot1=alotm+(1.-al)otp
[2788]	250	END IF
[3605]	251	ELSE
	252	IF (lfirst) ALLOCATE(v1(0,0,0,0))
[2788]	253	END IF
	254	!$OMP END MASTER
[3605]	255	!$OMP BARRIER
[2820]	256	IF(lfirst) THEN
[3086]	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)
[2820]	261	END IF
[3605]	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
[2788]	271	CALL scatter2d(v1,v2)
[3605]	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
[3086]	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(:)
[3605]	292
	293	!--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
	294	IF(.NOT.lAdjTro) THEN
[2788]	295	DO i=1,klon
[3086]	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(:)))
[2788]	301	END DO
[3086]	302	!-------------------------------------------------------------------------------
	303	ELSE !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
	304	!-------------------------------------------------------------------------------
[2788]	305	y_frac=(REAL(days_elapsed)+jH_cur)/year_len
[2968]	306
	307	!--- OUTPUT SIGMA LEVELS
[2788]	308	DO i=1,klon
	309
[3086]	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
[2788]	314
[3069]	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.
[2968]	318	SigT_in = get_SigTrop(i,itrp)
[3086]	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
[3069]	321
	322	!--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
[3086]	323	! dynamical tropopause (especially in filaments) are conterbalanced with
	324	! a filter ensuring it stays within a certain distance around input (file)
[3069]	325	! tropopause, hence avoiding avoid a too thick stretched region ; a final
	326	! extra-safety filter keeps the tropopause pressure value realistic.
[3086]	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)
[2968]	334
[3086]	335	!--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
[2968]	336	alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
	337
[3086]	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
[3141]	342	Sig_bot = MIN(Sig_bo0beta,0.1(9.+Sig_bo0)) !--- must be <1
[3086]	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
[2968]	352
[3086]	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]
[2788]	356	phi(:)=0.
[3086]	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)
[2788]	361
[3069]	362	!--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
[3086]	363	Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)*(alphaphi(:))
[2788]	364
[2968]	365	!--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
[3086]	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(:)))
[2968]	370
	371	!--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
[3086]	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', &
[2968]	374	5.E-9,3.0E-7)
	375	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
[3086]	376	ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
[2968]	377	5.E-8,1.5E-5)
	378	! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
	379
[3086]	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)
[3069]	388	END IF
[2788]	389	END DO
	390	END IF
[3086]	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
[3069]	402	CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
[3086]	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
[3069]	405	END IF
[2788]	406
	407	CONTAINS
	408
	409
	410	!-------------------------------------------------------------------------------
	411	!
	412	SUBROUTINE update_fields()
	413	!
	414	!-------------------------------------------------------------------------------
	415	IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION
[2968]	416	CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
[2820]	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.
[2788]	419	CALL get_3Dfields(v1) !--- Read ozone field(s)
[2820]	420	IF(lAdjTro) THEN !--- Additional files for fields strain
[3086]	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")
[2788]	424	END IF
	425	ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS
[2981]	426	IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
[2820]	427	IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
[2981]	428	CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
[2788]	429	IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS
[2981]	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)
[2788]	433	CALL get_3Dfields(v1m) !--- Read ozone field(s)
[2820]	434	IF(lAdjTro) THEN !--- Additional files for fields strain
[3086]	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")
[2788]	438	END IF
	439	ELSE !=== SHIFT FIELDS
[2981]	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)
[2788]	444	v1m=v1p !--- Ozone fields
[2820]	445	IF(lAdjTro) THEN !--- Additional files for fields strain
[3086]	446	IF(lPrSfile) pgm=pgp !--- Surface pressure
	447	IF(lPrTfile) ptm=ptp !--- Tropopause pressure
	448	IF(lO3Tfile) otm=otp !--- Tropopause ozone
[2788]	449	END IF
	450	END IF
	451	irec=irec+1
[2981]	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)
[2788]	455	CALL get_3Dfields(v1p) !--- Read ozone field(s)
[2820]	456	IF(lAdjTro) THEN !--- Additional files for fields strain
[3086]	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")
[2788]	460	END IF
[2981]	461	irec=irec-1
[2788]	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)
[2968]	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])
[2788]	486	CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
	487
	488	!--- Flip latitudes: ascending in input file, descending in "rlatu".
[2968]	489	IF(n_dim==3) THEN
[3086]	490	v(1,:) = v(1,nlat:1:-1)
	491	v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1) !--- Duplication
[2788]	492	ELSE
[3086]	493	v(:,:) = v(:,nlat:1:-1)
[2788]	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)
[2968]	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])
[2788]	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".
[2968]	522	IF(n_dim==3) THEN
[3086]	523	v(1,:,:,:) = v(1,nlat:1:-1,:,:)
	524	v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1) !--- Duplication
[2788]	525	ELSE
[3086]	526	v(:,:,:,:) = v(:,nlat:1:-1,:,:)
[2788]	527	END IF
	528
	529	END SUBROUTINE get_3Dfields
	530	!
	531	!-------------------------------------------------------------------------------
	532
	533
[2968]	534
[2788]	535	!-------------------------------------------------------------------------------
	536	!
[3086]	537	FUNCTION get_SigTrop(ih,it) RESULT(out)
[2788]	538	!
	539	!-------------------------------------------------------------------------------
	540	! Arguments:
[3086]	541	REAL :: out
[2788]	542	INTEGER, INTENT(IN) :: ih
	543	INTEGER, INTENT(OUT) :: it
	544	!-------------------------------------------------------------------------------
[3086]	545	!--- Pressure at tropopause read from the forcing file
	546	IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
[2968]	547
[3086]	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
[2968]	554	END FUNCTION get_SigTrop
	555	!
	556	!-------------------------------------------------------------------------------
	557
	558
	559	!-------------------------------------------------------------------------------
	560	!
[3086]	561	FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
[2968]	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:
[3086]	576	REAL :: out !--- Pressure at tropopause
[2968]	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
[3086]	580	REAL, INTENT(IN) :: pres(:) !--- Pressure profile, increasing
[2968]	581	REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV)
	582	REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause
	583	!-------------------------------------------------------------------------------
[2788]	584	! Local variables:
[3086]	585	REAL :: o3t !--- Ozone concent. at tropopause
	586	REAL :: al !--- Interpolation coefficient
	587	REAL :: coef !--- Coeff of latitude modulation
[2968]	588	REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause
[2788]	589	!-------------------------------------------------------------------------------
[2968]	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
[2788]	596	o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization
	597	END IF
	598
[2968]	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))
[3086]	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)
[2968]	605
[3086]	606	END FUNCTION PTrop_chem
[2788]	607	!
	608	!-------------------------------------------------------------------------------
	609
	610
[2968]	611	!-------------------------------------------------------------------------------
	612	!
[3086]	613	FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
[2968]	614	!
	615	!-------------------------------------------------------------------------------
	616	IMPLICIT NONE
	617	!-------------------------------------------------------------------------------
	618	! Arguments:
[3086]	619	LOGICAL :: out !--- .T. => some wrong values
[2968]	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
[3086]	635	out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
[2968]	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
[2788]	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