source: LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90 @ 4011

Last change on this file since 4011 was 2981, checked in by dcugnet, 7 years ago

Fix for time interpolation when ozone files have 14 records (monthly means).

File size: 29.6 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  IMPLICIT NONE
8
9!-------------------------------------------------------------------------------
10! Purpose: Regrid pressure with an averaging method. Operations done:
11!  * on the root process: read a NetCDF 3D or 4D field at a given day.
12!  * pack the fields to the LMDZ "horizontal "physics" grid and scatter
13!    to all threads of all processes;
14!  * in all the threads of all the processes, regrid the fields in pressure
15!    to the LMDZ vertical grid.
16!  * the forcing fields are stretched if the following arguments are present:
17!     - "lat_in":   input file latitudes.
18!     - "pcen_in":  input file cells center pressure levels.
19!     - "ptrop_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 boundaries.
37!  * Regridding in pressure is conservative, second order.
38!  * All the fields are regridded as a single multi-dimensional array, so it
39!    saves CPU time to call this procedure once for several NetCDF variables
40!    rather than several times, each time for a single NetCDF variable.
41!  * The input file pressure at tropopause can be (in decreasing priority):
42!    1) read from the file if "tropopause_air_pressure" field is available.
43!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
44!    3) computed using "tro3" and a fixed threshold otherwise, determined using
45!    an empirical three parameters law:
46!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
47!       => co1 and co2 are in ppbV, and co3 in degrees.
48!       => co3 allow a smooth transition between north and south hemispheres.
49!  * If available, the field "ps" (input file pressure at surface) is used.
50!    If not, the current LMDZ ground pressure is taken instead.
51!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
52!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
53!  * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
54!    interval "pint_in(k:k+1)]" and variable "nam(l)".
55!  * In the 2D file case, the values are the same for all longitudes.
56!  * The tropopause correction works like this: the input fields (file) are
57!  interpolated on output (LMDZ) pressure field, which is streched using a power
58!  law in a limited zone made of 3 layers:
59!    1) between the two tropopauses (file and LMDZ)
60!    2) in an upper and a lower transitions layers
61!  The stretching function has the following form:
62!    Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
63!  This value shifts the file tropopause to the height of the one of LMDZ.
64!  The stretching is fully applied in the central zone only, and only partially
65!  in the transitions zones, thick enough to guarantee a growing stretched
66!  pressure field. The ponderation function for alpha to modulate the stretching
67!  is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0)
68!  in the transition layers (from 1 to 0 ; in fact: sections of 1/log function),
69!  making this localisation function quasi-trapezoidal.
70!
71! * The following fields are on the global "dynamics" grid, as read from files:
72  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
73  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
74  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
75  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
76  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
77  INTEGER, SAVE :: ntim_in                         !--- Records nb in input file
78  INTEGER, SAVE :: itrp0                           !--- idx above chem tropop.
79  INTEGER, SAVE :: irec                            !--- Current time index
80!      * for daily   input files: current julian day number
81!      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
82  LOGICAL, SAVE :: linterp                         !--- Interpolation in time
83  LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
84  LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
85  LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
86  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
87!$OMP THREADPRIVATE(lfirst)
88  REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
89  REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
90  REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
91  REAL,    PARAMETER :: o3t0 = 1.E-7  !--- Nominal O3 vmr at tropopause
92  LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause
93
94CONTAINS
95
96!-------------------------------------------------------------------------------
97!
98SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
99                             time_in, lon_in, lat_in, pcen_in, ptrop_ou)
100!
101!-------------------------------------------------------------------------------
102  USE dimphy,         ONLY: klon
103  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
104                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
105  USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
106  USE assert_m,       ONLY: assert
107  USE assert_eq_m,    ONLY: assert_eq
108  USE comvert_mod,    ONLY: scaleheight
109  USE interpolation,  ONLY: locate
110  USE regr_conserv_m, ONLY: regr_conserv
111  USE regr_lint_m,    ONLY: regr_lint
112  USE slopes_m,       ONLY: slopes
113  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
114  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
115  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
116  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
117!-------------------------------------------------------------------------------
118! Arguments:
119  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
120  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
121  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
122  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
123  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
124  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
125  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
126                                              !    since Jan 1 of current year
127  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
128  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
129  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
130  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
131!-------------------------------------------------------------------------------
132! Local variables:
133  include "clesphys.h"
134  include "YOMCST.h"
135  CHARACTER(LEN=80)  :: sub
136  CHARACTER(LEN=320) :: str
137  INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
138  INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
139  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
140  REAL    :: y_frac                           !--- Elapsed year fraction
141  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
142  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
143  REAL    :: Sig_bot, Sig_top                 !--- Fully strained layer  bounds
144  REAL    :: Sig_bo0, Sig_to0                 !--- Total strained layers bounds
145  REAL    :: Sig_in(SIZE(pint_in))            !--- Input field sigma levels
146  REAL    :: Sig_ou (nbp_lev+1)               !--- Output LMDZ sigma levels
147  REAL    :: phi    (nbp_lev+1)               !--- Stretching exponent anomaly
148  REAL    :: pstr_ou(nbp_lev+1)               !--- Stretched pressure levels
149  REAL    :: pintou (nbp_lev+1)               !--- pint_ou in reversed order
150  REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
151  REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
152! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
153! v2: Field scattered to the partial "physics" horizontal grid.
154!     In the 2D file case, the values are the same for all longitudes.
155!     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
156! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
157  REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
158  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
159  LOGICAL :: ll
160!-------------------------------------------------------------------------------
161  sub="regr_pr_time_av"
162  nlev_in=SIZE(pint_in)-1
163  CALL assert(SIZE(v3,1)==klon,          TRIM(sub)//" v3 klon")
164  CALL assert(SIZE(v3,2)==nbp_lev,       TRIM(sub)//" v3 nbp_lev")
165  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
166  CALL assert(SIZE(pint_ou,1)==klon     ,TRIM(sub)//" pint_ou klon")
167  CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1")
168  IF(PRESENT(lon_in))   CALL assert(SIZE(lon_in  )==klon,TRIM(sub)//" lon_in klon")
169  IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
170  IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
171  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in")
172  lAdjTro=PRESENT(ptrop_ou)
173  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
174  CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
175
176  !$OMP MASTER
177  IF(is_mpi_root) THEN
178
179    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
180    IF(lfirst) THEN
181      lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
182      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
183      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
184      CALL NF95_INQ_DIMID(fID,"time",vID)
185      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
186      linterp=PRESENT(time_in).AND.ntim_in==14
187      IF(linterp) THEN
188        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
189        ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
190        ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
191        ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
192        IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
193      END IF
194      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
195      IF(lAdjTro) itrp0=locate(pcen_in,pTropUp)
196      CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub)
197      CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub)
198      CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub)
199      IF(lPrTrop) THEN
200        CALL msg(lAdjTro,'    PRESSURE AT TROPOPAUSE from file.')
201      ELSE IF(lO3Trop) THEN
202        CALL msg(lAdjTro,'    O3 CONCENTRATION AT TROPOPAUSE from file.')
203      ELSE IF(lo3tp) THEN
204        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
205      ELSE
206        CALL msg(lAdjTro,'    CONSTANT O3 concentration at tropopause.')
207      END IF
208    END IF
209
210    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
211    CALL update_fields()
212
213    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
214    IF(linterp) THEN
215      WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '&
216      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
217      CALL msg(.TRUE.,str,sub)
218      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
219      v1=al*v1m+(1.-al)*v1p
220      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
221      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
222      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
223    END IF
224  END IF
225  !$OMP END MASTER
226  IF(lfirst) THEN
227    lfirst=.FALSE.;      CALL bcast(lfirst)
228    IF(lAdjTro)          CALL bcast(itrp0)
229    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
230    CALL bcast(lO3Trop); CALL bcast(linterp)
231  END IF
232  CALL scatter2d(v1,v2)
233  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
234  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
235  IF(lPrTrop) CALL scatter2d(pt1,pt2)
236  IF(lO3Trop) CALL scatter2d(ot1,ot2)
237
238  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
239  IF(.NOT.lAdjTro) THEN
240    DO i=1,klon
241      pintou = pint_ou(i,nbp_lev+1:1:-1)
242      CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
243                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
244    END DO
245  ELSE
246    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
247
248    !--- OUTPUT SIGMA LEVELS
249    DO i=1,klon
250
251      !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
252      pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
253      Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
254      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
255
256      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
257      SigT_in = get_SigTrop(i,itrp)
258      SigT_ou = ptrop_ou(i)/ps2(i)
259
260      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
261      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
262      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
263      ptropou(i)=SigT_ou*ps2(i)
264
265      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
266      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
267
268      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
269      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
270      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
271
272      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
273      beta = LOG(Sig_top)/LOG(Sig_bot)
274      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
275      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
276
277      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
278      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
279      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
280      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
281      ibo0 = locate(Sig_ou(:),Sig_bo0)
282      ito0 = locate(Sig_ou(:),Sig_to0)
283
284      !--- FUNCTION FOR STRETCHING LOCALISATION
285      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
286      phi(:)=0.
287      phi(itop+1:ibot) =  1.
288      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
289                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
290      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
291                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
292
293      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
294      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
295
296      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
297      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
298                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
299
300      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
301      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
302      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
303                     5.E-9,3.0E-7)
304!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
305      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
306                     5.E-8,1.5E-5)
307!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
308
309    END DO
310  END IF
311
312CONTAINS
313
314
315!-------------------------------------------------------------------------------
316!
317SUBROUTINE update_fields()
318!
319!-------------------------------------------------------------------------------
320  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
321    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
322    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
323    !--- MIN -> Security in the unlikely case of roundup errors.
324    CALL get_3Dfields(v1)               !--- Read ozone field(s)
325    IF(lAdjTro) THEN                    !--- Additional files for fields strain
326      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
327      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
328      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
329    END IF
330  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
331    IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
332    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
333    CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
334    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
335      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&
336      //'reading record ',irec,' (time ',time_in(irec),')'
337      CALL msg(.TRUE.,str,sub)
338      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
339      IF(lAdjTro) THEN                  !--- Additional files for fields strain
340        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
341        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
342        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
343      END IF
344    ELSE                                !=== SHIFT FIELDS
345      irec=irec+1
346      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&
347      //' current next one (',irec,', time ',time_in(irec),')'
348      CALL msg(.TRUE.,str,sub)
349      v1m=v1p                           !--- Ozone fields
350      IF(lAdjTro) THEN                  !--- Additional files for fields strain
351        IF(lPrSurf) psm=psp             !--- Surface pressure
352        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
353        IF(lO3Trop) otm=otp             !--- Tropopause ozone
354      END IF
355    END IF
356    irec=irec+1
357    WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'&
358    ,irec,' (time ',time_in(irec),')'
359    CALL msg(.TRUE.,str,sub)
360    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
361    IF(lAdjTro) THEN                    !--- Additional files for fields strain
362      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
363      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
364      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
365    END IF
366    irec=irec-1
367  END IF
368
369END SUBROUTINE update_fields
370!
371!-------------------------------------------------------------------------------
372
373
374!-------------------------------------------------------------------------------
375!
376SUBROUTINE get_2Dfield(v,var)
377!
378!-------------------------------------------------------------------------------
379! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
380!          from the input file.
381! Remark: In case the field is zonal, it is duplicated along first dimension.
382!-------------------------------------------------------------------------------
383! Arguments:
384  REAL,             INTENT(INOUT) :: v(:,:)
385  CHARACTER(LEN=*), INTENT(IN)    :: var
386!-------------------------------------------------------------------------------
387  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
388  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
389  IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
390  IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
391  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
392
393  !--- Flip latitudes: ascending in input file, descending in "rlatu".
394  IF(n_dim==3) THEN
395    v(1,:) = v(1,nbp_lat:1:-1)
396    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
397  ELSE
398    v(:,:) = v(:,nbp_lat:1:-1)
399  END IF
400
401END SUBROUTINE get_2Dfield
402!
403!-------------------------------------------------------------------------------
404
405
406!-------------------------------------------------------------------------------
407!
408SUBROUTINE get_3Dfields(v)
409!
410!-------------------------------------------------------------------------------
411! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
412!          from the input file. Fields are stacked on fourth dimension.
413! Remark: In case the field is zonal, it is duplicated along first dimension.
414!-------------------------------------------------------------------------------
415! Arguments:
416  REAL, INTENT(INOUT) :: v(:,:,:,:)
417!-------------------------------------------------------------------------------
418  DO i=1,SIZE(nam)
419    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
420    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
421    IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
422    IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
423    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
424  END DO
425
426  !--- Flip latitudes: ascending in input file, descending in "rlatu".
427  IF(n_dim==3) THEN
428    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
429    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
430  ELSE
431    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
432  END IF
433
434END SUBROUTINE get_3Dfields
435!
436!-------------------------------------------------------------------------------
437
438
439
440!-------------------------------------------------------------------------------
441!
442FUNCTION get_SigTrop(ih,it)
443!
444!-------------------------------------------------------------------------------
445! Arguments:
446  INTEGER, INTENT(IN)  :: ih
447  INTEGER, INTENT(OUT) :: it
448  REAL                 :: get_Sigtrop
449!-------------------------------------------------------------------------------
450  !--- Pressure at tropopause is read in the forcing file
451  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
452    get_SigTrop=pt2(ih)/ps2(ih); RETURN
453  END IF
454  !--- Chemical tropopause definition is used using a particular threshold
455  IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
456    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
457  ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
458    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
459  ELSE                                         !--- o3trop CONSTANT
460    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
461  END IF
462  get_SigTrop=get_SigTrop/ps2(ih)
463
464END FUNCTION get_SigTrop
465!
466!-------------------------------------------------------------------------------
467
468
469!-------------------------------------------------------------------------------
470!
471FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
472!
473!-------------------------------------------------------------------------------
474! Purpose: Determine the tropopause using chemical criterium, ie the pressure
475!          at which the ozone concentration reaches a certain level.
476!-------------------------------------------------------------------------------
477! Remarks:
478! * Input field is upside down (increasing pressure // increasing vertical idx)
479!   The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
480!   where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
481!   the O3 profile in this zone is decreasing with pressure.
482! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
483!   determined using an empirical law roughly derived from ... & al.
484!-------------------------------------------------------------------------------
485! Arguments:
486  REAL ::                   chem_tropopause    !--- Pressure at tropopause
487  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
488  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
489  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
490  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
491  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
492  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
493  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
494!-------------------------------------------------------------------------------
495! Local variables:
496  REAL    :: o3t                               !--- Ozone concent. at tropopause
497  REAL    :: al                                !--- Interpolation coefficient
498  REAL    :: coef                              !--- Coeff of latitude modulation
499  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
500!-------------------------------------------------------------------------------
501  !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
502  IF(PRESENT(o3trop)) THEN                     !=== THRESHOLD FROM ARGUMENTS
503    o3t=o3trop
504  ELSE                                         !=== EMPIRICAL LAW
505    coef = TANH(lat_in(ih)/co3(3))             !--- Latitude  modulation
506    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Seasonal  modulation
507    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
508  END IF
509
510  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
511  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
512  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
513  IF(PRESENT(pcen)) THEN
514    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
515  ELSE
516    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
517  END IF
518  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
519
520END FUNCTION chem_tropopause
521!
522!-------------------------------------------------------------------------------
523
524
525!-------------------------------------------------------------------------------
526!
527FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
528!
529!-------------------------------------------------------------------------------
530  IMPLICIT NONE
531!-------------------------------------------------------------------------------
532! Arguments:
533  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
534  REAL,             INTENT(IN) :: o3col(:), lon, lat
535  INTEGER,          INTENT(IN) :: ilev0
536  CHARACTER(LEN=*), INTENT(IN) :: layer
537  REAL, OPTIONAL,   INTENT(IN) :: vmin
538  REAL, OPTIONAL,   INTENT(IN) :: vmax
539!-------------------------------------------------------------------------------
540! Local variables:
541  INTEGER :: k
542  LOGICAL :: lmin, lmax
543  REAL    :: fac
544  CHARACTER(LEN=6) :: unt
545!-------------------------------------------------------------------------------
546  !--- NOTHING TO DO
547  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
548  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
549  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
550
551  !--- SOME TOO LOW VALUES FOUND
552  IF(lmin) THEN
553    CALL unitc(vmin,unt,fac)
554    DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
555      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
556        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
557        fac*vmin,unt//' in '//TRIM(layer)
558    END DO
559  END IF
560
561  !--- SOME TOO HIGH VALUES FOUND
562  IF(lmax) THEN
563    CALL unitc(vmax,unt,fac)
564    DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
565      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
566        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
567        fac*vmax,unt//' in '//TRIM(layer)
568    END DO
569  END IF
570
571END FUNCTION check_ozone
572!
573!-------------------------------------------------------------------------------
574
575
576!-------------------------------------------------------------------------------
577!
578SUBROUTINE unitc(val,unt,fac)
579!
580!-------------------------------------------------------------------------------
581  IMPLICIT NONE
582!-------------------------------------------------------------------------------
583! Arguments:
584  REAL,             INTENT(IN)  :: val
585  CHARACTER(LEN=6), INTENT(OUT) :: unt
586  REAL,             INTENT(OUT) :: fac
587!-------------------------------------------------------------------------------
588! Local variables:
589  INTEGER :: ndgt
590!-------------------------------------------------------------------------------
591  ndgt=3*FLOOR(LOG10(val)/3.)
592  SELECT CASE(ndgt)
593    CASE(-6);     unt=' ppmV '; fac=1.E6
594    CASE(-9);     unt=' ppbV '; fac=1.E9
595    CASE DEFAULT; unt=' vmr  '; fac=1.0
596  END SELECT
597
598END SUBROUTINE unitc
599!
600!-------------------------------------------------------------------------------
601
602
603!-------------------------------------------------------------------------------
604!
605SUBROUTINE msg(ll,str,sub)
606!
607!-------------------------------------------------------------------------------
608  USE print_control_mod, ONLY: lunout
609  IMPLICIT NONE
610!-------------------------------------------------------------------------------
611! Arguments:
612  LOGICAL,                    INTENT(IN) :: ll
613  CHARACTER(LEN=*),           INTENT(IN) :: str
614  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
615!-------------------------------------------------------------------------------
616  IF(.NOT.ll) RETURN
617  IF(PRESENT(sub)) THEN
618    WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
619  ELSE
620    WRITE(lunout,*)TRIM(str)
621  END IF
622
623END SUBROUTINE msg
624!
625!-------------------------------------------------------------------------------
626
627END SUBROUTINE regr_pr_time_av
628!
629!-------------------------------------------------------------------------------
630
631END MODULE regr_pr_time_av_m
Note: See TracBrowser for help on using the repository browser.