source: LMDZ6/branches/Ocean_skin/libf/phylmd/regr_pr_time_av_m.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago

Merge revisions 3427:3600 of trunk into branch Ocean_skin

File size: 34.7 KB
RevLine 
[2788]1! $Id$
2MODULE regr_pr_time_av_m
3
4  USE write_field_phy
5  USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2968]6  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
[3069]7  USE print_control_mod,  ONLY: prt_level
[2788]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.
[2968]17!  * the forcing fields are stretched if the following arguments are present:
[3086]18!     - "lat_in":  input file latitudes.
19!     - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
[2968]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.
[2788]22!-------------------------------------------------------------------------------
23! Initial routine: regr_pr_av_m module (L. Guez).
24! Present version: David Cugnet ; corresponding additions:
25!    - time interpolation
26!    - 3D compliant
[2968]27!    - vertical stretching of the field to allow tropopause matching between
28!    input field and current lmdz field.
[2788]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.
[3086]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.
[2788]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.
[2968]48!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
[3086]49!    3) computed using "tro3" and a fixed threshold otherwise, constant or
50!    determined using an empirical three parameters law:
[2968]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.
[2788]54!  * If available, the field "ps" (input file pressure at surface) is used.
55!    If not, the current LMDZ ground pressure is taken instead.
[2968]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.
[3086]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)]".
[2968]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
[3086]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)
[2968]66!  The stretching function has the following form:
[3086]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.
[2968]75!
[2788]76! * The following fields are on the global "dynamics" grid, as read from files:
[3086]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
[2788]82  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
83  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
[2968]84  INTEGER, SAVE :: ntim_in                         !--- Records nb in input file
85  INTEGER, SAVE :: itrp0                           !--- idx above chem tropop.
[2788]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
[3086]90  LOGICAL, SAVE :: lPrSfile                        !--- Surface pressure flag
91  LOGICAL, SAVE :: lPrTfile                        !--- Tropopause pressure flag
92  LOGICAL, SAVE :: lO3Tfile                        !--- Tropopause ozone flag
[2788]93  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
[2819]94!$OMP THREADPRIVATE(lfirst)
[3069]95  REAL,    PARAMETER :: pTropUp=9.E+3 !--- Value  <  tropopause pressure (Pa)
[3086]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
[3087]100  LOGICAL, PARAMETER :: ldebug=.FALSE.!--- Force writefield_phy multiple outputs
[3069]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).
[2788]105
106CONTAINS
107
108!-------------------------------------------------------------------------------
109!
[3086]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)
[2788]112!
113!-------------------------------------------------------------------------------
114  USE dimphy,         ONLY: klon
[2968]115  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
116                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
[2788]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
[3605]120!!  USE comvert_mod,    ONLY: scaleheight
[2788]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
[3605]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
[2788]128  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
[3605]129  USE geometry_mod,                 ONLY: ind_cell_glo
[2788]130!-------------------------------------------------------------------------------
131! Arguments:
[3086]132  INTEGER,           INTENT(IN) :: fID        !--- NetCDF file ID
[2788]133  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
[3086]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)
[2788]142  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
143                                              !    since Jan 1 of current year
[3069]144  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
145  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
[3086]146  REAL,    INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres   (klon)
[2788]147!-------------------------------------------------------------------------------
148! Local variables:
149  include "clesphys.h"
150  include "YOMCST.h"
151  CHARACTER(LEN=80)  :: sub
[2968]152  CHARACTER(LEN=320) :: str
[3086]153  INTEGER :: vID, ncerr, n_var, ibot, iout, nn
154  INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
[2820]155  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
[2788]156  REAL    :: y_frac                           !--- Elapsed year fraction
[2968]157  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
158  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
[3086]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
[2788]171!     In the 2D file case, the values are the same for all longitudes.
[3086]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
[3605]178  INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
[2968]179  LOGICAL :: ll
[3069]180!--- For debug
[3086]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
[2788]184!-------------------------------------------------------------------------------
185  sub="regr_pr_time_av"
[3086]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")
[2968]193  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
[3086]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
[2788]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
[3086]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
[2968]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
[3086]218      ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
[2788]219      IF(linterp) THEN
[3086]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))
[2788]224      END IF
[2819]225      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
[3086]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')
[2788]235    END IF
236
237    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
[2820]238    CALL update_fields()
[2788]239
240    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
241    IF(linterp) THEN
[3086]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)
[2968]244      CALL msg(.TRUE.,str,sub)
[2788]245      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
246      v1=al*v1m+(1.-al)*v1p
[3086]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
[2788]250    END IF
[3605]251  ELSE
252    IF (lfirst) ALLOCATE(v1(0,0,0,0))
[2788]253  END IF
254  !$OMP END MASTER
[3605]255  !$OMP BARRIER
[2820]256  IF(lfirst) THEN
[3086]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)
[2820]261  END IF
[3605]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 
[2788]271  CALL scatter2d(v1,v2)
[3605]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
[3086]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(:)
[3605]292
293  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
294  IF(.NOT.lAdjTro) THEN
[2788]295    DO i=1,klon
[3086]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(:)))
[2788]301    END DO
[3086]302!-------------------------------------------------------------------------------
303  ELSE                        !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
304!-------------------------------------------------------------------------------
[2788]305    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
[2968]306
307    !--- OUTPUT SIGMA LEVELS
[2788]308    DO i=1,klon
309
[3086]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
[2788]314
[3069]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.
[2968]318      SigT_in = get_SigTrop(i,itrp)
[3086]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
[3069]321
322      !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
[3086]323      ! dynamical tropopause (especially in filaments) are conterbalanced with
324      ! a filter ensuring it stays within a certain distance around input (file)
[3069]325      ! tropopause, hence avoiding avoid a too thick stretched region ; a final
326      ! extra-safety filter keeps the tropopause pressure value realistic.
[3086]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)
[2968]334
[3086]335      !--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
[2968]336      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
337
[3086]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
[3141]342      Sig_bot = MIN(Sig_bo0*beta,0.1*(9.+Sig_bo0)) !--- must be <1
[3086]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
[2968]352
[3086]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]
[2788]356      phi(:)=0.
[3086]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)
[2788]361
[3069]362      !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
[3086]363      Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)**(alpha*phi(:))
[2788]364
[2968]365      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
[3086]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(:)))
[2968]370
371      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
[3086]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',  &
[2968]374                     5.E-9,3.0E-7)
375!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
[3086]376      ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
[2968]377                     5.E-8,1.5E-5)
378!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
379
[3086]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)
[3069]388      END IF
[2788]389    END DO
390  END IF
[3086]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
[3069]402    CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
[3086]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
[3069]405  END IF
[2788]406
407CONTAINS
408
409
410!-------------------------------------------------------------------------------
411!
412SUBROUTINE update_fields()
413!
414!-------------------------------------------------------------------------------
415  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
[2968]416    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
[2820]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.
[2788]419    CALL get_3Dfields(v1)               !--- Read ozone field(s)
[2820]420    IF(lAdjTro) THEN                    !--- Additional files for fields strain
[3086]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")
[2788]424    END IF
425  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
[2981]426    IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
[2820]427    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
[2981]428    CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
[2788]429    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
[2981]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)
[2788]433      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
[2820]434      IF(lAdjTro) THEN                  !--- Additional files for fields strain
[3086]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")
[2788]438      END IF
439    ELSE                                !=== SHIFT FIELDS
[2981]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)
[2788]444      v1m=v1p                           !--- Ozone fields
[2820]445      IF(lAdjTro) THEN                  !--- Additional files for fields strain
[3086]446        IF(lPrSfile) pgm=pgp             !--- Surface pressure
447        IF(lPrTfile) ptm=ptp             !--- Tropopause pressure
448        IF(lO3Tfile) otm=otp             !--- Tropopause ozone
[2788]449      END IF
450    END IF
451    irec=irec+1
[2981]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)
[2788]455    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
[2820]456    IF(lAdjTro) THEN                    !--- Additional files for fields strain
[3086]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")
[2788]460    END IF
[2981]461    irec=irec-1
[2788]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)
[2968]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])
[2788]486  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
487
488  !--- Flip latitudes: ascending in input file, descending in "rlatu".
[2968]489  IF(n_dim==3) THEN
[3086]490    v(1,:) = v(1,nlat:1:-1)
491    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1)  !--- Duplication
[2788]492  ELSE
[3086]493    v(:,:) = v(:,nlat:1:-1)
[2788]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)
[2968]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])
[2788]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".
[2968]522  IF(n_dim==3) THEN
[3086]523    v(1,:,:,:) = v(1,nlat:1:-1,:,:)
524    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1)  !--- Duplication
[2788]525  ELSE
[3086]526    v(:,:,:,:) = v(:,nlat:1:-1,:,:)
[2788]527  END IF
528
529END SUBROUTINE get_3Dfields
530!
531!-------------------------------------------------------------------------------
532
533
[2968]534
[2788]535!-------------------------------------------------------------------------------
536!
[3086]537FUNCTION get_SigTrop(ih,it) RESULT(out)
[2788]538!
539!-------------------------------------------------------------------------------
540! Arguments:
[3086]541  REAL                 :: out
[2788]542  INTEGER, INTENT(IN)  :: ih
543  INTEGER, INTENT(OUT) :: it
544!-------------------------------------------------------------------------------
[3086]545  !--- Pressure at tropopause read from the forcing file
546       IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
[2968]547
[3086]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
[2968]554END FUNCTION get_SigTrop
555!
556!-------------------------------------------------------------------------------
557
558
559!-------------------------------------------------------------------------------
560!
[3086]561FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
[2968]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:
[3086]576  REAL                        :: out           !--- Pressure at tropopause
[2968]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
[3086]580  REAL,           INTENT(IN)  :: pres(:)       !--- Pressure profile, increasing
[2968]581  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
582  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
583!-------------------------------------------------------------------------------
[2788]584! Local variables:
[3086]585  REAL :: o3t                                  !--- Ozone concent. at tropopause
586  REAL :: al                                   !--- Interpolation coefficient
587  REAL :: coef                                 !--- Coeff of latitude modulation
[2968]588  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
[2788]589!-------------------------------------------------------------------------------
[2968]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
[2788]596    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
597  END IF
598
[2968]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))
[3086]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)
[2968]605
[3086]606END FUNCTION PTrop_chem
[2788]607!
608!-------------------------------------------------------------------------------
609
610
[2968]611!-------------------------------------------------------------------------------
612!
[3086]613FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
[2968]614!
615!-------------------------------------------------------------------------------
616  IMPLICIT NONE
617!-------------------------------------------------------------------------------
618! Arguments:
[3086]619  LOGICAL                      :: out          !--- .T. => some wrong values
[2968]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
[3086]635  out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
[2968]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
[2788]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.