source: LMDZ6/branches/LMDZ-COSP/libf/phylmd/regr_pr_time_av_m.F90 @ 5409

Last change on this file since 5409 was 3598, checked in by jghattas, 5 years ago

Correction sur le commit [3596]. Desole.

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