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