[2788] | 1 | ! $Id$ |
---|
| 2 | MODULE 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 | |
---|
| 106 | CONTAINS |
---|
| 107 | |
---|
| 108 | !------------------------------------------------------------------------------- |
---|
| 109 | ! |
---|
[3086] | 110 | SUBROUTINE 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 |
---|
[4489] | 115 | USE netcdf95, ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, & |
---|
| 116 | NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, nf95_get_var |
---|
[5084] | 117 | USE netcdf, ONLY: NF90_INQ_VARID, NF90_NOERR |
---|
[2788] | 118 | USE assert_m, ONLY: assert |
---|
| 119 | USE assert_eq_m, ONLY: assert_eq |
---|
[3435] | 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 |
---|
[3435] | 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 |
---|
[3435] | 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 |
---|
[3435] | 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 |
---|
[3596] | 251 | ELSE |
---|
[3598] | 252 | IF (lfirst) ALLOCATE(v1(0,0,0,0)) |
---|
[2788] | 253 | END IF |
---|
| 254 | !$OMP END MASTER |
---|
[3479] | 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 |
---|
[3435] | 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) |
---|
[3435] | 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(:) |
---|
[3435] | 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 | |
---|
| 407 | CONTAINS |
---|
| 408 | |
---|
| 409 | |
---|
| 410 | !------------------------------------------------------------------------------- |
---|
| 411 | ! |
---|
| 412 | SUBROUTINE 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 | |
---|
| 464 | END SUBROUTINE update_fields |
---|
| 465 | ! |
---|
| 466 | !------------------------------------------------------------------------------- |
---|
| 467 | |
---|
| 468 | |
---|
| 469 | !------------------------------------------------------------------------------- |
---|
| 470 | ! |
---|
| 471 | SUBROUTINE 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) |
---|
[4489] | 484 | IF(n_dim==2) call NF95_GET_VAR(fID,vID,v(1,:), start=[ 1,irec]) |
---|
| 485 | IF(n_dim==3) call NF95_GET_VAR(fID,vID,v(:,:), start=[1,1,irec]) |
---|
[2788] | 486 | |
---|
| 487 | !--- Flip latitudes: ascending in input file, descending in "rlatu". |
---|
[2968] | 488 | IF(n_dim==3) THEN |
---|
[3086] | 489 | v(1,:) = v(1,nlat:1:-1) |
---|
| 490 | v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1) !--- Duplication |
---|
[2788] | 491 | ELSE |
---|
[3086] | 492 | v(:,:) = v(:,nlat:1:-1) |
---|
[2788] | 493 | END IF |
---|
| 494 | |
---|
| 495 | END SUBROUTINE get_2Dfield |
---|
| 496 | ! |
---|
| 497 | !------------------------------------------------------------------------------- |
---|
| 498 | |
---|
| 499 | |
---|
| 500 | !------------------------------------------------------------------------------- |
---|
| 501 | ! |
---|
| 502 | SUBROUTINE get_3Dfields(v) |
---|
| 503 | ! |
---|
| 504 | !------------------------------------------------------------------------------- |
---|
| 505 | ! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam" |
---|
| 506 | ! from the input file. Fields are stacked on fourth dimension. |
---|
| 507 | ! Remark: In case the field is zonal, it is duplicated along first dimension. |
---|
| 508 | !------------------------------------------------------------------------------- |
---|
| 509 | ! Arguments: |
---|
| 510 | REAL, INTENT(INOUT) :: v(:,:,:,:) |
---|
| 511 | !------------------------------------------------------------------------------- |
---|
| 512 | DO i=1,SIZE(nam) |
---|
| 513 | CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID) |
---|
[2968] | 514 | CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim) |
---|
[4489] | 515 | IF(n_dim==3) call NF95_GET_VAR(fID,vID,v(1,:,:,i), start=[ 1,1,irec]) |
---|
| 516 | IF(n_dim==4) call NF95_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec]) |
---|
[2788] | 517 | END DO |
---|
| 518 | |
---|
| 519 | !--- Flip latitudes: ascending in input file, descending in "rlatu". |
---|
[2968] | 520 | IF(n_dim==3) THEN |
---|
[3086] | 521 | v(1,:,:,:) = v(1,nlat:1:-1,:,:) |
---|
| 522 | v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1) !--- Duplication |
---|
[2788] | 523 | ELSE |
---|
[3086] | 524 | v(:,:,:,:) = v(:,nlat:1:-1,:,:) |
---|
[2788] | 525 | END IF |
---|
| 526 | |
---|
| 527 | END SUBROUTINE get_3Dfields |
---|
| 528 | ! |
---|
| 529 | !------------------------------------------------------------------------------- |
---|
| 530 | |
---|
| 531 | |
---|
[2968] | 532 | |
---|
[2788] | 533 | !------------------------------------------------------------------------------- |
---|
| 534 | ! |
---|
[3086] | 535 | FUNCTION get_SigTrop(ih,it) RESULT(out) |
---|
[2788] | 536 | ! |
---|
| 537 | !------------------------------------------------------------------------------- |
---|
| 538 | ! Arguments: |
---|
[3086] | 539 | REAL :: out |
---|
[2788] | 540 | INTEGER, INTENT(IN) :: ih |
---|
| 541 | INTEGER, INTENT(OUT) :: it |
---|
| 542 | !------------------------------------------------------------------------------- |
---|
[3086] | 543 | !--- Pressure at tropopause read from the forcing file |
---|
| 544 | IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF |
---|
[2968] | 545 | |
---|
[3086] | 546 | !--- Chemical tropopause definition based on a particular threshold |
---|
| 547 | IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih)) |
---|
| 548 | ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1)) |
---|
| 549 | ELSE ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF |
---|
| 550 | out=out/Pgnd_in(ih) |
---|
| 551 | |
---|
[2968] | 552 | END FUNCTION get_SigTrop |
---|
| 553 | ! |
---|
| 554 | !------------------------------------------------------------------------------- |
---|
| 555 | |
---|
| 556 | |
---|
| 557 | !------------------------------------------------------------------------------- |
---|
| 558 | ! |
---|
[3086] | 559 | FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out) |
---|
[2968] | 560 | ! |
---|
| 561 | !------------------------------------------------------------------------------- |
---|
| 562 | ! Purpose: Determine the tropopause using chemical criterium, ie the pressure |
---|
| 563 | ! at which the ozone concentration reaches a certain level. |
---|
| 564 | !------------------------------------------------------------------------------- |
---|
| 565 | ! Remarks: |
---|
| 566 | ! * Input field is upside down (increasing pressure // increasing vertical idx) |
---|
| 567 | ! The sweep is done from top to ground, starting at the 50hPa layer (idx it0), |
---|
| 568 | ! where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached: |
---|
| 569 | ! the O3 profile in this zone is decreasing with pressure. |
---|
| 570 | ! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or |
---|
| 571 | ! determined using an empirical law roughly derived from ... & al. |
---|
| 572 | !------------------------------------------------------------------------------- |
---|
| 573 | ! Arguments: |
---|
[3086] | 574 | REAL :: out !--- Pressure at tropopause |
---|
[2968] | 575 | INTEGER, INTENT(IN) :: ih !--- Horizontal index |
---|
| 576 | INTEGER, INTENT(OUT) :: it !--- Index of tropopause layer |
---|
| 577 | INTEGER, INTENT(IN) :: it0 !--- Idx: higher than tropopause |
---|
[3086] | 578 | REAL, INTENT(IN) :: pres(:) !--- Pressure profile, increasing |
---|
[2968] | 579 | REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV) |
---|
| 580 | REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause |
---|
| 581 | !------------------------------------------------------------------------------- |
---|
[2788] | 582 | ! Local variables: |
---|
[3086] | 583 | REAL :: o3t !--- Ozone concent. at tropopause |
---|
| 584 | REAL :: al !--- Interpolation coefficient |
---|
| 585 | REAL :: coef !--- Coeff of latitude modulation |
---|
[2968] | 586 | REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause |
---|
[2788] | 587 | !------------------------------------------------------------------------------- |
---|
[2968] | 588 | !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN |
---|
| 589 | IF(PRESENT(o3trop)) THEN !=== THRESHOLD FROM ARGUMENTS |
---|
| 590 | o3t=o3trop |
---|
| 591 | ELSE !=== EMPIRICAL LAW |
---|
| 592 | coef = TANH(lat_in(ih)/co3(3)) !--- Latitude modulation |
---|
| 593 | coef = SIN (2.*RPI*(y_frac-1./6.)) * coef !--- Seasonal modulation |
---|
[2788] | 594 | o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization |
---|
| 595 | END IF |
---|
| 596 | |
---|
[2968] | 597 | !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1) |
---|
| 598 | it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO |
---|
| 599 | al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1)) |
---|
[3086] | 600 | IF(Ploc=='C') out = pres(it)**(1.-al) * pres(it+1)**al |
---|
| 601 | IF(Ploc=='I') out = SQRT(pres(it)**(1.-al) * pres(it+2)**al *pres(it+1)) |
---|
| 602 | it = locate(pres(:), out) !--- pres(it)<Ptrop<pres(it+1) |
---|
[2968] | 603 | |
---|
[3086] | 604 | END FUNCTION PTrop_chem |
---|
[2788] | 605 | ! |
---|
| 606 | !------------------------------------------------------------------------------- |
---|
| 607 | |
---|
| 608 | |
---|
[2968] | 609 | !------------------------------------------------------------------------------- |
---|
| 610 | ! |
---|
[3086] | 611 | FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out) |
---|
[2968] | 612 | ! |
---|
| 613 | !------------------------------------------------------------------------------- |
---|
| 614 | IMPLICIT NONE |
---|
| 615 | !------------------------------------------------------------------------------- |
---|
| 616 | ! Arguments: |
---|
[3086] | 617 | LOGICAL :: out !--- .T. => some wrong values |
---|
[2968] | 618 | REAL, INTENT(IN) :: o3col(:), lon, lat |
---|
| 619 | INTEGER, INTENT(IN) :: ilev0 |
---|
| 620 | CHARACTER(LEN=*), INTENT(IN) :: layer |
---|
| 621 | REAL, OPTIONAL, INTENT(IN) :: vmin |
---|
| 622 | REAL, OPTIONAL, INTENT(IN) :: vmax |
---|
| 623 | !------------------------------------------------------------------------------- |
---|
| 624 | ! Local variables: |
---|
| 625 | INTEGER :: k |
---|
| 626 | LOGICAL :: lmin, lmax |
---|
| 627 | REAL :: fac |
---|
| 628 | CHARACTER(LEN=6) :: unt |
---|
| 629 | !------------------------------------------------------------------------------- |
---|
| 630 | !--- NOTHING TO DO |
---|
| 631 | lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0 |
---|
| 632 | lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0 |
---|
[3086] | 633 | out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN |
---|
[2968] | 634 | |
---|
| 635 | !--- SOME TOO LOW VALUES FOUND |
---|
| 636 | IF(lmin) THEN |
---|
| 637 | CALL unitc(vmin,unt,fac) |
---|
| 638 | DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE |
---|
| 639 | WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' & |
---|
| 640 | //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', & |
---|
| 641 | fac*vmin,unt//' in '//TRIM(layer) |
---|
| 642 | END DO |
---|
| 643 | END IF |
---|
| 644 | |
---|
| 645 | !--- SOME TOO HIGH VALUES FOUND |
---|
| 646 | IF(lmax) THEN |
---|
| 647 | CALL unitc(vmax,unt,fac) |
---|
| 648 | DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE |
---|
| 649 | WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' & |
---|
| 650 | //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', & |
---|
| 651 | fac*vmax,unt//' in '//TRIM(layer) |
---|
| 652 | END DO |
---|
| 653 | END IF |
---|
| 654 | |
---|
| 655 | END FUNCTION check_ozone |
---|
| 656 | ! |
---|
| 657 | !------------------------------------------------------------------------------- |
---|
| 658 | |
---|
| 659 | |
---|
| 660 | !------------------------------------------------------------------------------- |
---|
| 661 | ! |
---|
| 662 | SUBROUTINE unitc(val,unt,fac) |
---|
| 663 | ! |
---|
| 664 | !------------------------------------------------------------------------------- |
---|
| 665 | IMPLICIT NONE |
---|
| 666 | !------------------------------------------------------------------------------- |
---|
| 667 | ! Arguments: |
---|
| 668 | REAL, INTENT(IN) :: val |
---|
| 669 | CHARACTER(LEN=6), INTENT(OUT) :: unt |
---|
| 670 | REAL, INTENT(OUT) :: fac |
---|
| 671 | !------------------------------------------------------------------------------- |
---|
| 672 | ! Local variables: |
---|
| 673 | INTEGER :: ndgt |
---|
| 674 | !------------------------------------------------------------------------------- |
---|
| 675 | ndgt=3*FLOOR(LOG10(val)/3.) |
---|
| 676 | SELECT CASE(ndgt) |
---|
| 677 | CASE(-6); unt=' ppmV '; fac=1.E6 |
---|
| 678 | CASE(-9); unt=' ppbV '; fac=1.E9 |
---|
| 679 | CASE DEFAULT; unt=' vmr '; fac=1.0 |
---|
| 680 | END SELECT |
---|
| 681 | |
---|
| 682 | END SUBROUTINE unitc |
---|
| 683 | ! |
---|
| 684 | !------------------------------------------------------------------------------- |
---|
| 685 | |
---|
| 686 | |
---|
| 687 | !------------------------------------------------------------------------------- |
---|
| 688 | ! |
---|
| 689 | SUBROUTINE msg(ll,str,sub) |
---|
| 690 | ! |
---|
| 691 | !------------------------------------------------------------------------------- |
---|
| 692 | USE print_control_mod, ONLY: lunout |
---|
| 693 | IMPLICIT NONE |
---|
| 694 | !------------------------------------------------------------------------------- |
---|
| 695 | ! Arguments: |
---|
| 696 | LOGICAL, INTENT(IN) :: ll |
---|
| 697 | CHARACTER(LEN=*), INTENT(IN) :: str |
---|
| 698 | CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub |
---|
| 699 | !------------------------------------------------------------------------------- |
---|
| 700 | IF(.NOT.ll) RETURN |
---|
| 701 | IF(PRESENT(sub)) THEN |
---|
| 702 | WRITE(lunout,*)TRIM(sub)//': '//TRIM(str) |
---|
| 703 | ELSE |
---|
| 704 | WRITE(lunout,*)TRIM(str) |
---|
| 705 | END IF |
---|
| 706 | |
---|
| 707 | END SUBROUTINE msg |
---|
| 708 | ! |
---|
| 709 | !------------------------------------------------------------------------------- |
---|
| 710 | |
---|
[2788] | 711 | END SUBROUTINE regr_pr_time_av |
---|
| 712 | ! |
---|
| 713 | !------------------------------------------------------------------------------- |
---|
| 714 | |
---|
| 715 | END MODULE regr_pr_time_av_m |
---|