source: LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/regr_pr_time_av_m.F90 @ 2980

Last change on this file since 2980 was 2980, checked in by acozic, 7 years ago

merges with trunk for rev 2968 and 2971 --
2968 = improved method for ozone forcing
2971 = control outputs for debug are removed

File size: 29.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  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 paramztrized 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
204        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
205      END IF
206    END IF
207
208    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
209    CALL update_fields()
210
211    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
212    IF(linterp) THEN
213      WRITE(str,'(3(a,f7.3))')'Interpolating O3 at julian day ',julien,' from '&
214      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
215      CALL msg(.TRUE.,str,sub)
216      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
217      v1=al*v1m+(1.-al)*v1p
218      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
219      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
220      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
221    END IF
222  END IF
223  !$OMP END MASTER
224  IF(lfirst) THEN
225    lfirst=.FALSE.;      CALL bcast(lfirst)
226    IF(lAdjTro)          CALL bcast(itrp0)
227    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
228    CALL bcast(lO3Trop); CALL bcast(linterp)
229  END IF
230  CALL scatter2d(v1,v2)
231  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
232  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
233  IF(lPrTrop) CALL scatter2d(pt1,pt2)
234  IF(lO3Trop) CALL scatter2d(ot1,ot2)
235
236  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
237  IF(.NOT.lAdjTro) THEN
238    DO i=1,klon
239      pintou = pint_ou(i,nbp_lev+1:1:-1)
240      CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
241                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
242    END DO
243  ELSE
244    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
245
246    !--- OUTPUT SIGMA LEVELS
247    DO i=1,klon
248
249      !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
250      pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
251      Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
252      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
253
254      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
255      SigT_in = get_SigTrop(i,itrp)
256      SigT_ou = ptrop_ou(i)/ps2(i)
257
258      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
259      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
260      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
261      ptropou(i)=SigT_ou*ps2(i)
262
263      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
264      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
265
266      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
267      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
268      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
269
270      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
271      beta = LOG(Sig_top)/LOG(Sig_bot)
272      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
273      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
274
275      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
276      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
277      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
278      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
279      ibo0 = locate(Sig_ou(:),Sig_bo0)
280      ito0 = locate(Sig_ou(:),Sig_to0)
281
282      !--- FUNCTION FOR STRETCHING LOCALISATION
283      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
284      phi(:)=0.
285      phi(itop+1:ibot) =  1.
286      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
287                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
288      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
289                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
290
291      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
292      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
293
294      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
295      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
296                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
297
298      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
299      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
300      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
301                     5.E-9,3.0E-7)
302!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
303      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
304                     5.E-8,1.5E-5)
305!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
306
307    END DO
308  END IF
309
310CONTAINS
311
312
313!-------------------------------------------------------------------------------
314!
315SUBROUTINE update_fields()
316!
317!-------------------------------------------------------------------------------
318  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
319    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
320    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
321    !--- MIN -> Security in the unlikely case of roundup errors.
322    CALL get_3Dfields(v1)               !--- Read ozone field(s)
323    IF(lAdjTro) THEN                    !--- Additional files for fields strain
324      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
325      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
326      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
327    END IF
328  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
329    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
330    CALL msg(.TRUE.,sub,'Refreshing adjacent Ozone forcing fields.')
331    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
332      irec=locate(time_in,julien)       !--- Need to locate surrounding times
333      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
334      IF(lAdjTro) THEN                  !--- Additional files for fields strain
335        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
336        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
337        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
338      END IF
339    ELSE                                !=== SHIFT FIELDS
340      v1m=v1p                           !--- Ozone fields
341      IF(lAdjTro) THEN                  !--- Additional files for fields strain
342        IF(lPrSurf) psm=psp             !--- Surface pressure
343        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
344        IF(lO3Trop) otm=otp             !--- Tropopause ozone
345      END IF
346    END IF
347    irec=irec+1
348    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
349    IF(lAdjTro) THEN                    !--- Additional files for fields strain
350      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
351      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
352      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
353    END IF
354    IF(lfirst) irec=irec-1
355  END IF
356
357END SUBROUTINE update_fields
358!
359!-------------------------------------------------------------------------------
360
361
362!-------------------------------------------------------------------------------
363!
364SUBROUTINE get_2Dfield(v,var)
365!
366!-------------------------------------------------------------------------------
367! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
368!          from the input file.
369! Remark: In case the field is zonal, it is duplicated along first dimension.
370!-------------------------------------------------------------------------------
371! Arguments:
372  REAL,             INTENT(INOUT) :: v(:,:)
373  CHARACTER(LEN=*), INTENT(IN)    :: var
374!-------------------------------------------------------------------------------
375  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
376  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
377  IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
378  IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
379  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
380
381  !--- Flip latitudes: ascending in input file, descending in "rlatu".
382  IF(n_dim==3) THEN
383    v(1,:) = v(1,nbp_lat:1:-1)
384    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
385  ELSE
386    v(:,:) = v(:,nbp_lat:1:-1)
387  END IF
388
389END SUBROUTINE get_2Dfield
390!
391!-------------------------------------------------------------------------------
392
393
394!-------------------------------------------------------------------------------
395!
396SUBROUTINE get_3Dfields(v)
397!
398!-------------------------------------------------------------------------------
399! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
400!          from the input file. Fields are stacked on fourth dimension.
401! Remark: In case the field is zonal, it is duplicated along first dimension.
402!-------------------------------------------------------------------------------
403! Arguments:
404  REAL, INTENT(INOUT) :: v(:,:,:,:)
405!-------------------------------------------------------------------------------
406  DO i=1,SIZE(nam)
407    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
408    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
409    IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
410    IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
411    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
412  END DO
413
414  !--- Flip latitudes: ascending in input file, descending in "rlatu".
415  IF(n_dim==3) THEN
416    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
417    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
418  ELSE
419    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
420  END IF
421
422END SUBROUTINE get_3Dfields
423!
424!-------------------------------------------------------------------------------
425
426
427
428!-------------------------------------------------------------------------------
429!
430FUNCTION get_SigTrop(ih,it)
431!
432!-------------------------------------------------------------------------------
433! Arguments:
434  INTEGER, INTENT(IN)  :: ih
435  INTEGER, INTENT(OUT) :: it
436  REAL                 :: get_Sigtrop
437!-------------------------------------------------------------------------------
438  !--- Pressure at tropopause is read in the forcing file
439  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
440    get_SigTrop=pt2(ih)/ps2(ih); RETURN
441  END IF
442  !--- Chemical tropopause definition is used using a particular threshold
443  IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
444    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
445  ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
446    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
447  ELSE                                         !--- o3trop CONSTANT
448    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
449  END IF
450  get_SigTrop=get_SigTrop/ps2(ih)
451
452END FUNCTION get_SigTrop
453!
454!-------------------------------------------------------------------------------
455
456
457!-------------------------------------------------------------------------------
458!
459FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
460!
461!-------------------------------------------------------------------------------
462! Purpose: Determine the tropopause using chemical criterium, ie the pressure
463!          at which the ozone concentration reaches a certain level.
464!-------------------------------------------------------------------------------
465! Remarks:
466! * Input field is upside down (increasing pressure // increasing vertical idx)
467!   The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
468!   where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
469!   the O3 profile in this zone is decreasing with pressure.
470! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
471!   determined using an empirical law roughly derived from ... & al.
472!-------------------------------------------------------------------------------
473! Arguments:
474  REAL ::                   chem_tropopause    !--- Pressure at tropopause
475  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
476  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
477  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
478  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
479  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
480  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
481  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
482!-------------------------------------------------------------------------------
483! Local variables:
484  REAL    :: o3t                               !--- Ozone concent. at tropopause
485  REAL    :: al                                !--- Interpolation coefficient
486  REAL    :: coef                              !--- Coeff of latitude modulation
487  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
488!-------------------------------------------------------------------------------
489  !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
490  IF(PRESENT(o3trop)) THEN                     !=== THRESHOLD FROM ARGUMENTS
491    o3t=o3trop
492  ELSE                                         !=== EMPIRICAL LAW
493    coef = TANH(lat_in(ih)/co3(3))             !--- Latitude  modulation
494    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Seasonal  modulation
495    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
496  END IF
497
498  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
499  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
500  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
501  IF(PRESENT(pcen)) THEN
502    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
503  ELSE
504    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
505  END IF
506  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
507
508END FUNCTION chem_tropopause
509!
510!-------------------------------------------------------------------------------
511
512
513!-------------------------------------------------------------------------------
514!
515FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
516!
517!-------------------------------------------------------------------------------
518  IMPLICIT NONE
519!-------------------------------------------------------------------------------
520! Arguments:
521  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
522  REAL,             INTENT(IN) :: o3col(:), lon, lat
523  INTEGER,          INTENT(IN) :: ilev0
524  CHARACTER(LEN=*), INTENT(IN) :: layer
525  REAL, OPTIONAL,   INTENT(IN) :: vmin
526  REAL, OPTIONAL,   INTENT(IN) :: vmax
527!-------------------------------------------------------------------------------
528! Local variables:
529  INTEGER :: k
530  LOGICAL :: lmin, lmax
531  REAL    :: fac
532  CHARACTER(LEN=6) :: unt
533!-------------------------------------------------------------------------------
534  !--- NOTHING TO DO
535  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
536  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
537  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
538
539  !--- SOME TOO LOW VALUES FOUND
540  IF(lmin) THEN
541    CALL unitc(vmin,unt,fac)
542    DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
543      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
544        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
545        fac*vmin,unt//' in '//TRIM(layer)
546    END DO
547  END IF
548
549  !--- SOME TOO HIGH VALUES FOUND
550  IF(lmax) THEN
551    CALL unitc(vmax,unt,fac)
552    DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
553      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
554        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
555        fac*vmax,unt//' in '//TRIM(layer)
556    END DO
557  END IF
558
559END FUNCTION check_ozone
560!
561!-------------------------------------------------------------------------------
562
563
564!-------------------------------------------------------------------------------
565!
566SUBROUTINE unitc(val,unt,fac)
567!
568!-------------------------------------------------------------------------------
569  IMPLICIT NONE
570!-------------------------------------------------------------------------------
571! Arguments:
572  REAL,             INTENT(IN)  :: val
573  CHARACTER(LEN=6), INTENT(OUT) :: unt
574  REAL,             INTENT(OUT) :: fac
575!-------------------------------------------------------------------------------
576! Local variables:
577  INTEGER :: ndgt
578!-------------------------------------------------------------------------------
579  ndgt=3*FLOOR(LOG10(val)/3.)
580  SELECT CASE(ndgt)
581    CASE(-6);     unt=' ppmV '; fac=1.E6
582    CASE(-9);     unt=' ppbV '; fac=1.E9
583    CASE DEFAULT; unt=' vmr  '; fac=1.0
584  END SELECT
585
586END SUBROUTINE unitc
587!
588!-------------------------------------------------------------------------------
589
590
591!-------------------------------------------------------------------------------
592!
593SUBROUTINE msg(ll,str,sub)
594!
595!-------------------------------------------------------------------------------
596  USE print_control_mod, ONLY: lunout
597  IMPLICIT NONE
598!-------------------------------------------------------------------------------
599! Arguments:
600  LOGICAL,                    INTENT(IN) :: ll
601  CHARACTER(LEN=*),           INTENT(IN) :: str
602  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
603!-------------------------------------------------------------------------------
604  IF(.NOT.ll) RETURN
605  IF(PRESENT(sub)) THEN
606    WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
607  ELSE
608    WRITE(lunout,*)TRIM(str)
609  END IF
610
611END SUBROUTINE msg
612!
613!-------------------------------------------------------------------------------
614
615END SUBROUTINE regr_pr_time_av
616!
617!-------------------------------------------------------------------------------
618
619END MODULE regr_pr_time_av_m
Note: See TracBrowser for help on using the repository browser.