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

Last change on this file since 2966 was 2820, checked in by dcugnet, 8 years ago

Fix and cleaning/improvements over the previous commit:

  • Old style climoz files are handled without trouble.
  • New style climoz (monthly with 14 records or daily) are handled

with the same read_climoz key values (1 or 2).

  • Daily files reading activation is done using the same key

used in ce0l to create them: ok_daily_climoz=y. This key is now
put to FALSE by default (old style 360 days files are not considered
as daily).

File size: 19.7 KB
Line 
1! $Id$
2MODULE regr_pr_time_av_m
3
4  USE print_control_mod, ONLY: lunout
5  USE write_field_phy
6  USE mod_phys_lmdz_transfert_para, ONLY: bcast
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!  * input files fields are stretched in the viscinity (+/- 5 kms) of the mean
17!    tropopause (geometrical mean of LMDZ and input fields tropopauses) to force
18!    the resulting field tropopause height to the one of LMDZ. To switch this
19!    feature on, the following arguments must be present:
20!     - "lat_in":   input file latitudes.
21!     - "pcen_in":  input file cells center pressure levels.
22!     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
23!    Note that the ozone quantity conservation is not ensured.
24!-------------------------------------------------------------------------------
25! Initial routine: regr_pr_av_m module (L. Guez).
26! Present version: David Cugnet ; corresponding additions:
27!    - time interpolation
28!    - 3D compliant
29!    - vertical stretching of the field to allow tropopause and ground pressure
30!      matching between input field and current lmdz field.
31!-------------------------------------------------------------------------------
32! Remarks:
33!  * 3D fields are zonal means, with dimensions (latitude, pressure, julian day)
34!  * 4D fields have the dimensions:  (longitude, latitude, pressure, julian day)
35!  * All the fields are already on the horizontal grid (rlatu) or (rlatu,rlonv),
36!    except that the latitudes are in ascending order in the input file.
37!  * We assume that all the input fields have the same coordinates.
38!  * The target vertical LMDZ grid is the grid of layer boundaries.
39!  * Regridding in pressure is conservative, second order.
40!  * All the fields are regridded as a single multi-dimensional array, so it
41!    saves CPU time to call this procedure once for several NetCDF variables
42!    rather than several times, each time for a single NetCDF variable.
43!  * The input file pressure at tropopause can be (in decreasing priority):
44!    1) read from the file if "tropopause_air_pressure" field is available.
45!    2) computed from the input file ozone field using:
46!       - o3 concentration at tropopause if "tro3_at_tropopause" is available.
47!       - a default value (100ppbv) if not.
48!  * If available, the field "ps" (input file pressure at surface) is used.
49!    If not, the current LMDZ ground pressure is taken instead.
50!  * The O3 threshold for tropopause is defined using 3 coefficients:
51!    o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/c)
52!     => co1 and co2 are in ppmV, and co3 in degrees.
53!     => co3 allow a smooth transition between north and south hemispheres.
54!-------------------------------------------------------------------------------
55! * Fields with suffix "m"/"p" are at the closest records earlier/later than the
56!   middle of the julian day "julien", on the global "dynamics" horizontal grid.
57! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
58!   interval "pint_in(k:k+1)]" and variable "nam(l)".
59! * In the 2D file case, the values are the same for all longitudes.
60! * The following fields are on the global "dynamics" grid, as read from files:
61  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
62  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
63  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
64  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
65  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
66  INTEGER, SAVE :: irec                            !--- Current time index
67!      * for daily   input files: current julian day number
68!      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
69  LOGICAL, SAVE :: linterp                         !--- Interpolation in time
70  LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
71  LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
72  LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
73  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
74!$OMP THREADPRIVATE(lfirst)
75  REAL, PARAMETER :: co3(3)=[91.,28.,20.]          !--- Coeffs for o3 threshold
76  REAL, PARAMETER :: prtrop=5.E+3  !--- Value lower than the tropopause pressure
77  REAL, PARAMETER :: delta =5.     !--- Dist. around tropopause for strain (kms)
78
79CONTAINS
80
81!-------------------------------------------------------------------------------
82!
83SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, time_in,    &
84                                     lat_in, pcen_in, ptrop_ou)
85!
86!-------------------------------------------------------------------------------
87  USE dimphy,         ONLY: klon
88  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err
89  USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
90  USE assert_m,       ONLY: assert
91  USE assert_eq_m,    ONLY: assert_eq
92  USE comvert_mod,    ONLY: scaleheight
93  USE interpolation,  ONLY: locate
94  USE regr_conserv_m, ONLY: regr_conserv
95  USE regr_lint_m,    ONLY: regr_lint
96  USE slopes_m,       ONLY: slopes
97  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
98  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
99  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
100  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
101!-------------------------------------------------------------------------------
102! Arguments:
103  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
104  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
105  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
106  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
107  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ Pres, Pa (klon,llm+1)
108  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
109  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
110                                              !    since Jan 1 of current year
111  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
112  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
113  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pressure (klon)
114!-------------------------------------------------------------------------------
115! Local variables:
116  include "clesphys.h"
117  include "YOMCST.h"
118  CHARACTER(LEN=80)  :: sub
119  CHARACTER(LEN=320) :: msg
120  INTEGER :: vID, ncerr, n_var, nlev_in,ntim_in, ndim, i, ibot, itop, itrp,itrp0
121  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
122  REAL    :: y_frac                           !--- Elapsed year fraction
123  REAL    :: alpha, beta, al                  !--- For strectching/interpolation
124  REAL    :: SigT_in, SigT_ou, SigT_mn        !--- Tropopause: in/out/mean
125  REAL    :: SigA_bot, SigA_top               !--- Strained domain bounds
126  REAL    :: Sig_in (SIZE(pint_in))           !--- Input field sigma levels
127  REAL    :: phi    (SIZE(pint_in))           !--- Stretching exponent anomaly
128  REAl    :: Pint_st(SIZE(pint_in))           !--- Stretched pressure levels
129  REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
130  REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
131! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
132! v2: Field scattered to the partial "physics" horizontal grid.
133!     In the 2D file case, the values are the same for all longitudes.
134!     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
135! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
136  REAL, DIMENSION(nbp_lon, nbp_lat) ::  ps1, pt1, ot1
137  REAL, DIMENSION(klon)             ::  ps2, pt2, ot2
138!-------------------------------------------------------------------------------
139  sub="regr_pr_time_av"
140  nlev_in=SIZE(pint_in)-1; ntim_in=SIZE(time_in)
141  CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
142  CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
143  n_var = assert_eq(SIZE(nam), SIZE(v3,3),      TRIM(sub)//" v3 n_var")
144  CALL assert(SHAPE(pint_ou)==[klon, nbp_lev+1],TRIM(sub)//" pint_ou")
145  IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
146  IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
147  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in)==nlev_in,TRIM(sub)//" pcen_in")
148  lAdjTro=PRESENT(ptrop_ou)
149  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
150  CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
151
152  !$OMP MASTER
153  IF(is_mpi_root) THEN
154
155    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
156    IF(lfirst) THEN
157      lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
158      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
159      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
160      linterp=PRESENT(time_in); IF(linterp) linterp=ntim_in==14
161      IF(linterp) THEN
162        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
163        ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
164        ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
165        ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
166        IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
167      END IF
168      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
169      IF(lAdjTro) itrp0=locate(pcen_in,prtrop)
170      IF(lPrSurf) WRITE(lunout,*)TRIM(sub)//': Using GROUND PRESSURE from input O3 forcing file.'
171      IF(linterp) WRITE(lunout,*)TRIM(sub)//': Monthly O3 files => ONLINE TIME INTERPOLATION.'
172      IF(lAdjTro) THEN
173        WRITE(lunout,*)TRIM(sub)//': o3 forcing file tropopause location uses:'
174        IF(lPrTrop) THEN
175          WRITE(lunout,*)'    PRESSURE AT TROPOPAUSE from file.'
176        ELSE IF(lO3Trop) THEN
177          WRITE(lunout,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
178        ELSE
179          WRITE(lunout,*)'    PARAMETRIZED O3 concentration at tropopause.'
180        END IF
181      END IF
182    END IF
183
184    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
185    CALL update_fields()
186
187
188    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
189    IF(linterp) THEN
190      WRITE(lunout,'(3(a,f7.3))')TRIM(sub)//': Interpolating O3 at julian day '&
191      &,julien,' from fields at times ',time_in(irec),' and ', time_in(irec+1)
192      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
193      v1=al*v1m+(1.-al)*v1p
194      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
195      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
196      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
197    END IF
198  END IF
199  !$OMP END MASTER
200  IF(lfirst) THEN
201    lfirst=.FALSE. ;     CALL bcast(lfirst) ; CALL bcast(lPrSurf)
202    CALL bcast(lPrTrop); CALL bcast(lO3Trop); CALL bcast(linterp)
203    IF(lAdjTro) CALL bcast(itrp0)
204  END IF
205  CALL scatter2d(v1,v2)
206  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
207  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=Pint_ou(:,1); END IF
208  IF(lPrTrop) CALL scatter2d(pt1,pt2)
209  IF(lO3Trop) CALL scatter2d(ot1,ot2)
210
211  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
212  IF(.NOT.lAdjTro) THEN
213    DO i=1,klon
214      CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1),    &
215                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in))
216    END DO
217  ELSE
218    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
219    DO i=1,klon
220      SigT_in = get_SigTrop(i,itrp)        !--- input  (file) tropopause
221      SigT_ou = ptrop_ou(i)/pint_ou(i,1)   !--- output (lmdz) tropopause
222      SigT_mn = SQRT(SigT_in*SigT_ou)      !--- mean tropopause>strained domain
223
224      !--- DEFINE THE VERTICAL LAYER IN WHICH THE PROFILE IS STRECHED
225      beta=EXP(delta/scaleheight); Sig_in(:) = [pint_in(1:nlev_in)/ps2(i),1.]
226      SigA_bot = SigT_mn*beta ; ibot=locate(Sig_in(:),SigA_bot)
227      SigA_top = SigT_mn/beta ; itop=locate(Sig_in(:),SigA_top)+1
228
229      !--- HAT FUNCTION phi (/=0 in [SigA_bot,SigA_top] only)
230      phi(:)=0.
231      phi(itop:itrp)=(Sig_in(itop:itrp)-SigA_top)/(SigT_in-SigA_top)
232      phi(itrp:ibot)=(SigA_bot-Sig_in(itrp:ibot))/(SigA_bot-SigT_in)
233
234      !--- STRAINED INPUT logP PROFILE ; alpha: MAX. STRETCH. EXPONENT INCREMENT
235      alpha = LOG(SigT_ou/SigT_in)/LOG(SigT_in)
236      Pint_st(:) = pint_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
237
238      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL LEVELS
239      CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pint_ou(i,nbp_lev+1:1:-1),     &
240                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),Pint_st))
241    END DO
242  END IF
243
244
245CONTAINS
246
247
248!-------------------------------------------------------------------------------
249!
250SUBROUTINE update_fields()
251!
252!-------------------------------------------------------------------------------
253  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
254    WRITE(lunout,*)TRIM(sub)//': Updating Ozone forcing field: read from file.'
255    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
256    !--- MIN -> Security in the unlikely case of roundup errors.
257    CALL get_3Dfields(v1)               !--- Read ozone field(s)
258    IF(lAdjTro) THEN                    !--- Additional files for fields strain
259      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
260      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
261      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
262    END IF
263  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
264    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
265    WRITE(lunout,*)TRIM(sub)//': Refreshing adjacent Ozone forcing fields.'
266    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
267      irec=locate(time_in,julien)       !--- Need to locate surrounding times
268      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
269      IF(lAdjTro) THEN                  !--- Additional files for fields strain
270        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
271        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
272        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
273      END IF
274    ELSE                                !=== SHIFT FIELDS
275      v1m=v1p                           !--- Ozone fields
276      IF(lAdjTro) THEN                  !--- Additional files for fields strain
277        IF(lPrSurf) psm=psp             !--- Surface pressure
278        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
279        IF(lO3Trop) otm=otp             !--- Tropopause ozone
280      END IF
281    END IF
282    irec=irec+1
283    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
284    IF(lAdjTro) THEN                    !--- Additional files for fields strain
285      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
286      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
287      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
288    END IF
289    IF(lfirst) irec=irec-1
290  END IF
291
292END SUBROUTINE update_fields
293!
294!-------------------------------------------------------------------------------
295
296
297!-------------------------------------------------------------------------------
298!
299SUBROUTINE get_2Dfield(v,var)
300!
301!-------------------------------------------------------------------------------
302! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
303!          from the input file.
304! Remark: In case the field is zonal, it is duplicated along first dimension.
305!-------------------------------------------------------------------------------
306! Arguments:
307  REAL,             INTENT(INOUT) :: v(:,:)
308  CHARACTER(LEN=*), INTENT(IN)    :: var
309!-------------------------------------------------------------------------------
310  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
311  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
312  IF(ndim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
313  IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
314  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
315
316  !--- Flip latitudes: ascending in input file, descending in "rlatu".
317  IF(ndim==3) THEN
318    v(1,:) = v(1,nbp_lat:1:-1)
319    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
320  ELSE
321    v(:,:) = v(:,nbp_lat:1:-1)
322  END IF
323
324END SUBROUTINE get_2Dfield
325!
326!-------------------------------------------------------------------------------
327
328
329!-------------------------------------------------------------------------------
330!
331SUBROUTINE get_3Dfields(v)
332!
333!-------------------------------------------------------------------------------
334! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
335!          from the input file. Fields are stacked on fourth dimension.
336! Remark: In case the field is zonal, it is duplicated along first dimension.
337!-------------------------------------------------------------------------------
338! Arguments:
339  REAL, INTENT(INOUT) :: v(:,:,:,:)
340!-------------------------------------------------------------------------------
341  DO i=1,SIZE(nam)
342    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
343    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
344    IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
345    IF(ndim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
346    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
347  END DO
348
349  !--- Flip latitudes: ascending in input file, descending in "rlatu".
350  IF(ndim==3) THEN
351    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
352    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
353  ELSE
354    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
355  END IF
356
357END SUBROUTINE get_3Dfields
358!
359!-------------------------------------------------------------------------------
360
361
362!-------------------------------------------------------------------------------
363!
364FUNCTION get_SigTrop(ih,it)
365!
366!-------------------------------------------------------------------------------
367! Arguments:
368  INTEGER, INTENT(IN)  :: ih
369  INTEGER, INTENT(OUT) :: it
370  REAL                 :: get_Sigtrop
371!-------------------------------------------------------------------------------
372! Local variables:
373  INTEGER :: ii                                !--- Idx of layer containing o3t
374  REAL    :: o3t                               !--- Ozone concent. at tropopause
375  REAL    :: prt                               !--- Air pressure   at tropopause
376  REAL    :: al                                !--- Interpolation coefficient
377  REAL    :: coef                              !--- Coef: North/South transition
378!-------------------------------------------------------------------------------
379  !--- DETERMINE PRESSURE AT TROPOPAUSE AND DIVIDE IT BY GROUND PRESSURE
380  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
381    get_SigTrop=pt2(ih)/ps2(ih)
382    it=locate(pint_in(:),pt2(ih))
383  ELSE                                         !--- O3 THRESHOLD
384    coef = TANH(lat_in(i)/co3(3))              !--- Latitude dependant ponderat.
385    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Time-dependant ponderation
386    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
387    IF(lO3Trop) o3t=ot2(ih)                    !--- Value from file
388    !--- Starts from 50hPa and go down.
389    it=itrp0; DO WHILE(v2(ih,it+1,1)>=o3t); it=it+1; END DO
390    al=(v2(ih,it,1)-o3t)/(v2(ih,it,1)-v2(ih,it+1,1))
391    get_SigTrop = ( pcen_in(it)**(1.-al) * pcen_in(it+1)**al )/ps2(ih)
392  END IF
393
394END FUNCTION get_SigTrop
395!
396!-------------------------------------------------------------------------------
397
398
399END SUBROUTINE regr_pr_time_av
400!
401!-------------------------------------------------------------------------------
402
403END MODULE regr_pr_time_av_m
Note: See TracBrowser for help on using the repository browser.