source: LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90 @ 3141

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