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