[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 |
---|
[2788] | 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. |
---|
[2968] | 16 | ! * the forcing fields are stretched if the following arguments are present: |
---|
[3411] | 17 | ! - "lat_in": input file latitudes. |
---|
| 18 | ! - "pcen_in": input file cells center pressure levels. |
---|
| 19 | ! - "ptrop_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. |
---|
[3411] | 36 | ! * The target vertical LMDZ grid is the grid of layer boundaries. |
---|
| 37 | ! * Regridding in pressure is conservative, second order. |
---|
[2788] | 38 | ! * All the fields are regridded as a single multi-dimensional array, so it |
---|
| 39 | ! saves CPU time to call this procedure once for several NetCDF variables |
---|
| 40 | ! rather than several times, each time for a single NetCDF variable. |
---|
| 41 | ! * The input file pressure at tropopause can be (in decreasing priority): |
---|
| 42 | ! 1) read from the file if "tropopause_air_pressure" field is available. |
---|
[2968] | 43 | ! 2) computed using "tro3" and "tro3_at_tropopause' (if available). |
---|
[3411] | 44 | ! 3) computed using "tro3" and a fixed threshold otherwise, determined using |
---|
| 45 | ! an empirical three parameters law: |
---|
[2968] | 46 | ! o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3) |
---|
| 47 | ! => co1 and co2 are in ppbV, and co3 in degrees. |
---|
| 48 | ! => co3 allow a smooth transition between north and south hemispheres. |
---|
[2788] | 49 | ! * If available, the field "ps" (input file pressure at surface) is used. |
---|
| 50 | ! If not, the current LMDZ ground pressure is taken instead. |
---|
[2968] | 51 | ! * Fields with suffix "m"/"p" are at the closest records earlier/later than |
---|
| 52 | ! the mid-julian day "julien", on the global "dynamics" horizontal grid. |
---|
[3411] | 53 | ! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure |
---|
| 54 | ! interval "pint_in(k:k+1)]" and variable "nam(l)". |
---|
[2968] | 55 | ! * In the 2D file case, the values are the same for all longitudes. |
---|
| 56 | ! * The tropopause correction works like this: the input fields (file) are |
---|
| 57 | ! interpolated on output (LMDZ) pressure field, which is streched using a power |
---|
[3411] | 58 | ! law in a limited zone made of 3 layers: |
---|
| 59 | ! 1) between the two tropopauses (file and LMDZ) |
---|
| 60 | ! 2) in an upper and a lower transitions layers |
---|
[2968] | 61 | ! The stretching function has the following form: |
---|
[3411] | 62 | ! Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou) |
---|
| 63 | ! This value shifts the file tropopause to the height of the one of LMDZ. |
---|
| 64 | ! The stretching is fully applied in the central zone only, and only partially |
---|
| 65 | ! in the transitions zones, thick enough to guarantee a growing stretched |
---|
| 66 | ! pressure field. The ponderation function for alpha to modulate the stretching |
---|
| 67 | ! is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0) |
---|
| 68 | ! in the transition layers (from 1 to 0 ; in fact: sections of 1/log function), |
---|
| 69 | ! making this localisation function quasi-trapezoidal. |
---|
[2968] | 70 | ! |
---|
[2788] | 71 | ! * The following fields are on the global "dynamics" grid, as read from files: |
---|
[3411] | 72 | REAL, SAVE, ALLOCATABLE :: v1m(:,:,:,:) !--- Previous ozone fields |
---|
| 73 | REAL, SAVE, ALLOCATABLE :: v1p(:,:,:,:) !--- Next ozone fields |
---|
| 74 | REAL, SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure |
---|
[2788] | 75 | REAL, SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure |
---|
| 76 | REAL, SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio |
---|
[2968] | 77 | INTEGER, SAVE :: ntim_in !--- Records nb in input file |
---|
| 78 | INTEGER, SAVE :: itrp0 !--- idx above chem tropop. |
---|
[2788] | 79 | INTEGER, SAVE :: irec !--- Current time index |
---|
| 80 | ! * for daily input files: current julian day number |
---|
| 81 | ! * for monthly input files: julien is in [time_in(irec),time_in(irec+1)] |
---|
| 82 | LOGICAL, SAVE :: linterp !--- Interpolation in time |
---|
[3411] | 83 | LOGICAL, SAVE :: lPrSurf !--- Surface pressure flag |
---|
| 84 | LOGICAL, SAVE :: lPrTrop !--- Tropopause pressure flag |
---|
| 85 | LOGICAL, SAVE :: lO3Trop !--- Tropopause ozone flag |
---|
[2788] | 86 | LOGICAL, SAVE :: lfirst=.TRUE. !--- First call flag |
---|
[2819] | 87 | !$OMP THREADPRIVATE(lfirst) |
---|
[3411] | 88 | REAL, PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa) |
---|
| 89 | REAL, PARAMETER :: gamm = 0.4 !--- Relative thickness of transitions |
---|
| 90 | REAL, PARAMETER :: rho = 1.3 !--- Max tropopauses sigma ratio |
---|
| 91 | REAL, PARAMETER :: o3t0 = 1.E-7 !--- Nominal O3 vmr at tropopause |
---|
| 92 | LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause |
---|
[2788] | 93 | |
---|
| 94 | CONTAINS |
---|
| 95 | |
---|
| 96 | !------------------------------------------------------------------------------- |
---|
| 97 | ! |
---|
[3411] | 98 | SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, & |
---|
| 99 | time_in, lon_in, lat_in, pcen_in, ptrop_ou) |
---|
[2788] | 100 | ! |
---|
| 101 | !------------------------------------------------------------------------------- |
---|
| 102 | USE dimphy, ONLY: klon |
---|
[2968] | 103 | USE netcdf95, ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, & |
---|
| 104 | NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION |
---|
[2788] | 105 | USE netcdf, ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR |
---|
| 106 | USE assert_m, ONLY: assert |
---|
| 107 | USE assert_eq_m, ONLY: assert_eq |
---|
[3413] | 108 | !! USE comvert_mod, ONLY: scaleheight |
---|
[2788] | 109 | USE interpolation, ONLY: locate |
---|
| 110 | USE regr_conserv_m, ONLY: regr_conserv |
---|
| 111 | USE regr_lint_m, ONLY: regr_lint |
---|
| 112 | USE slopes_m, ONLY: slopes |
---|
[3413] | 113 | USE mod_phys_lmdz_para, ONLY: is_mpi_root,is_master |
---|
| 114 | USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid_type, unstructured |
---|
| 115 | USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter, gather |
---|
[2788] | 116 | USE phys_cal_mod, ONLY: calend, year_len, days_elapsed, jH_cur |
---|
[3413] | 117 | USE geometry_mod, ONLY: ind_cell_glo |
---|
[2788] | 118 | !------------------------------------------------------------------------------- |
---|
| 119 | ! Arguments: |
---|
[3411] | 120 | INTEGER, INTENT(IN) :: fID !--- NetCDF file ID |
---|
[2788] | 121 | CHARACTER(LEN=13), INTENT(IN) :: nam(:) !--- NetCDF variables names |
---|
[3411] | 122 | REAL, INTENT(IN) :: julien !--- Days since Jan 1st |
---|
| 123 | REAL, INTENT(IN) :: pint_in(:) !--- Interfaces file Pres, Pa, ascending |
---|
| 124 | REAL, INTENT(IN) :: pint_ou(:,:) !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1) |
---|
| 125 | REAL, INTENT(OUT) :: v3(:,:,:) !--- Regridded field (klon,llm,SIZE(nam)) |
---|
[2788] | 126 | REAL, INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days |
---|
| 127 | ! since Jan 1 of current year |
---|
[3411] | 128 | REAL, INTENT(IN), OPTIONAL :: lon_in(:) !--- Input/output longitudes vector |
---|
| 129 | REAL, INTENT(IN), OPTIONAL :: lat_in(:) !--- Input/output latitudes vector |
---|
| 130 | REAL, INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure |
---|
| 131 | REAL, INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon) |
---|
[2788] | 132 | !------------------------------------------------------------------------------- |
---|
| 133 | ! Local variables: |
---|
| 134 | include "clesphys.h" |
---|
| 135 | include "YOMCST.h" |
---|
| 136 | CHARACTER(LEN=80) :: sub |
---|
[2968] | 137 | CHARACTER(LEN=320) :: str |
---|
[3411] | 138 | INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp |
---|
| 139 | INTEGER :: i, nlev_in, n_dim, itop, ito0, i0 |
---|
[2820] | 140 | LOGICAL :: lAdjTro !--- Need to adjust tropopause |
---|
[2788] | 141 | REAL :: y_frac !--- Elapsed year fraction |
---|
[2968] | 142 | REAL :: alpha, beta, al !--- For stretching/interpolation |
---|
| 143 | REAL :: SigT_in, SigT_ou !--- Input and output tropopauses |
---|
[3411] | 144 | REAL :: Sig_bot, Sig_top !--- Fully strained layer bounds |
---|
| 145 | REAL :: Sig_bo0, Sig_to0 !--- Total strained layers bounds |
---|
| 146 | REAL :: Sig_in(SIZE(pint_in)) !--- Input field sigma levels |
---|
| 147 | REAL :: Sig_ou (nbp_lev+1) !--- Output LMDZ sigma levels |
---|
| 148 | REAL :: phi (nbp_lev+1) !--- Stretching exponent anomaly |
---|
| 149 | REAL :: pstr_ou(nbp_lev+1) !--- Stretched pressure levels |
---|
| 150 | REAL :: pintou (nbp_lev+1) !--- pint_ou in reversed order |
---|
| 151 | REAL :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam)) |
---|
| 152 | REAL :: v2(klon, SIZE(pint_in)-1,SIZE(nam)) |
---|
| 153 | ! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid. |
---|
| 154 | ! v2: Field scattered to the partial "physics" horizontal grid. |
---|
[2788] | 155 | ! In the 2D file case, the values are the same for all longitudes. |
---|
[3411] | 156 | ! "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)". |
---|
| 157 | ! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)" |
---|
| 158 | REAL, DIMENSION(nbp_lon, nbp_lat) :: ps1, pt1, ot1 |
---|
| 159 | REAL, DIMENSION(klon) :: ps2, pt2, ot2, ptropou |
---|
[3413] | 160 | INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:) |
---|
[2968] | 161 | LOGICAL :: ll |
---|
[2788] | 162 | !------------------------------------------------------------------------------- |
---|
| 163 | sub="regr_pr_time_av" |
---|
[3411] | 164 | nlev_in=SIZE(pint_in)-1 |
---|
| 165 | CALL assert(SIZE(v3,1)==klon, TRIM(sub)//" v3 klon") |
---|
| 166 | CALL assert(SIZE(v3,2)==nbp_lev, TRIM(sub)//" v3 nbp_lev") |
---|
[2968] | 167 | n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var") |
---|
[3411] | 168 | CALL assert(SIZE(pint_ou,1)==klon ,TRIM(sub)//" pint_ou klon") |
---|
| 169 | CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1") |
---|
| 170 | IF(PRESENT(lon_in)) CALL assert(SIZE(lon_in )==klon,TRIM(sub)//" lon_in klon") |
---|
| 171 | IF(PRESENT(lat_in)) CALL assert(SIZE(lat_in )==klon,TRIM(sub)//" lat_in klon") |
---|
| 172 | IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon") |
---|
| 173 | IF(PRESENT(pcen_in)) CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in") |
---|
| 174 | lAdjTro=PRESENT(ptrop_ou) |
---|
| 175 | IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) & |
---|
| 176 | CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1) |
---|
[2788] | 177 | |
---|
| 178 | !$OMP MASTER |
---|
| 179 | IF(is_mpi_root) THEN |
---|
| 180 | |
---|
| 181 | !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE |
---|
| 182 | IF(lfirst) THEN |
---|
[3411] | 183 | lPrSurf=NF90_INQ_VARID(fID,"ps" ,vID)==NF90_NOERR |
---|
| 184 | lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR |
---|
| 185 | lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause" ,vID)==NF90_NOERR |
---|
[2968] | 186 | CALL NF95_INQ_DIMID(fID,"time",vID) |
---|
| 187 | CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in) |
---|
| 188 | linterp=PRESENT(time_in).AND.ntim_in==14 |
---|
[2788] | 189 | IF(linterp) THEN |
---|
[3411] | 190 | ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var)) |
---|
| 191 | ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var)) |
---|
| 192 | ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat)) |
---|
| 193 | ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat)) |
---|
| 194 | IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat)) |
---|
[2788] | 195 | END IF |
---|
[2819] | 196 | !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa) |
---|
[3411] | 197 | IF(lAdjTro) itrp0=locate(pcen_in,pTropUp) |
---|
| 198 | CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub) |
---|
| 199 | CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub) |
---|
| 200 | CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub) |
---|
| 201 | IF(lPrTrop) THEN |
---|
| 202 | CALL msg(lAdjTro,' PRESSURE AT TROPOPAUSE from file.') |
---|
| 203 | ELSE IF(lO3Trop) THEN |
---|
| 204 | CALL msg(lAdjTro,' O3 CONCENTRATION AT TROPOPAUSE from file.') |
---|
| 205 | ELSE IF(lo3tp) THEN |
---|
| 206 | CALL msg(lAdjTro,' PARAMETRIZED O3 concentration at tropopause.') |
---|
| 207 | ELSE |
---|
| 208 | CALL msg(lAdjTro,' CONSTANT O3 concentration at tropopause.') |
---|
| 209 | END IF |
---|
[2788] | 210 | END IF |
---|
| 211 | |
---|
| 212 | !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES) |
---|
[2820] | 213 | CALL update_fields() |
---|
[2788] | 214 | |
---|
| 215 | !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES |
---|
| 216 | IF(linterp) THEN |
---|
[3411] | 217 | WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '& |
---|
| 218 | &//'fields at times ',time_in(irec),' and ', time_in(irec+1) |
---|
[2968] | 219 | CALL msg(.TRUE.,str,sub) |
---|
[2788] | 220 | al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec)) |
---|
| 221 | v1=al*v1m+(1.-al)*v1p |
---|
[3411] | 222 | IF(lPrSurf) ps1=al*psm+(1.-al)*psp |
---|
| 223 | IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp |
---|
| 224 | IF(lO3Trop) ot1=al*otm+(1.-al)*otp |
---|
[2788] | 225 | END IF |
---|
| 226 | END IF |
---|
| 227 | !$OMP END MASTER |
---|
[2820] | 228 | IF(lfirst) THEN |
---|
[3411] | 229 | lfirst=.FALSE.; CALL bcast(lfirst) |
---|
| 230 | IF(lAdjTro) CALL bcast(itrp0) |
---|
| 231 | CALL bcast(lPrTrop); CALL bcast(lPrSurf) |
---|
| 232 | CALL bcast(lO3Trop); CALL bcast(linterp) |
---|
[2820] | 233 | END IF |
---|
[3413] | 234 | |
---|
| 235 | IF (is_master) THEN |
---|
| 236 | ALLOCATE(ind_cell_glo_glo(klon_glo)) |
---|
| 237 | ELSE |
---|
| 238 | ALLOCATE(ind_cell_glo_glo(0)) |
---|
| 239 | ENDIF |
---|
| 240 | CALL gather(ind_cell_glo,ind_cell_glo_glo) |
---|
| 241 | IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:) |
---|
| 242 | |
---|
[2788] | 243 | CALL scatter2d(v1,v2) |
---|
[3413] | 244 | |
---|
[3411] | 245 | !--- No "ps" in input file => assumed to be equal to current LMDZ ground press |
---|
[3413] | 246 | IF(lPrSurf) THEN |
---|
| 247 | IF (is_master .AND. grid_type==unstructured) ps1(:,:)=ps1(:,ind_cell_glo_glo(:)) |
---|
| 248 | CALL scatter2d(ps1,ps2) |
---|
| 249 | ELSE |
---|
| 250 | ps2=pint_ou(:,1) |
---|
| 251 | END IF |
---|
[3411] | 252 | |
---|
[3413] | 253 | IF(lPrTrop) THEN |
---|
| 254 | IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:)) |
---|
| 255 | CALL scatter2d(pt1,pt2) |
---|
| 256 | ENDIF |
---|
| 257 | |
---|
| 258 | IF(lO3Trop) THEN |
---|
| 259 | IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:)) |
---|
| 260 | CALL scatter2d(ot1,ot2) |
---|
| 261 | ENDIF |
---|
[3411] | 262 | !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing |
---|
| 263 | IF(.NOT.lAdjTro) THEN |
---|
[2788] | 264 | DO i=1,klon |
---|
[3411] | 265 | pintou = pint_ou(i,nbp_lev+1:1:-1) |
---|
| 266 | CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:), & |
---|
| 267 | v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:))) |
---|
[2788] | 268 | END DO |
---|
[3411] | 269 | ELSE |
---|
[2788] | 270 | y_frac=(REAL(days_elapsed)+jH_cur)/year_len |
---|
[2968] | 271 | |
---|
| 272 | !--- OUTPUT SIGMA LEVELS |
---|
[2788] | 273 | DO i=1,klon |
---|
| 274 | |
---|
[3411] | 275 | !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES |
---|
| 276 | pintou=pint_ou(i,nbp_lev+1:1:-1) !--- increasing values |
---|
| 277 | Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)] !--- increasing values |
---|
| 278 | Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values |
---|
[2788] | 279 | |
---|
[3411] | 280 | !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE |
---|
[2968] | 281 | SigT_in = get_SigTrop(i,itrp) |
---|
[3411] | 282 | SigT_ou = ptrop_ou(i)/ps2(i) |
---|
[2968] | 283 | |
---|
[3411] | 284 | !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION |
---|
| 285 | IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho |
---|
| 286 | IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho |
---|
| 287 | ptropou(i)=SigT_ou*ps2(i) |
---|
[2968] | 288 | |
---|
[3411] | 289 | !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW |
---|
[2968] | 290 | alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou) |
---|
| 291 | |
---|
[3411] | 292 | !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES) |
---|
| 293 | Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot) |
---|
| 294 | Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top) |
---|
[2968] | 295 | |
---|
[3411] | 296 | !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE |
---|
| 297 | beta = LOG(Sig_top)/LOG(Sig_bot) |
---|
| 298 | Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta) |
---|
| 299 | Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top ** beta |
---|
| 300 | |
---|
| 301 | !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS |
---|
| 302 | !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1 |
---|
| 303 | Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0) |
---|
| 304 | Sig_to0 = Sig_to0*EXP(-gamm*ABS(alpha)) |
---|
| 305 | ibo0 = locate(Sig_ou(:),Sig_bo0) |
---|
| 306 | ito0 = locate(Sig_ou(:),Sig_to0) |
---|
| 307 | |
---|
| 308 | !--- FUNCTION FOR STRETCHING LOCALISATION |
---|
| 309 | ! 0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1 |
---|
[2788] | 310 | phi(:)=0. |
---|
[3411] | 311 | phi(itop+1:ibot) = 1. |
---|
| 312 | phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))& |
---|
| 313 | *LOG(Sig_top)/LOG(Sig_top/Sig_to0) |
---|
| 314 | phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))& |
---|
| 315 | *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0) |
---|
[2788] | 316 | |
---|
[3411] | 317 | !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER) |
---|
| 318 | pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:)) |
---|
[2788] | 319 | |
---|
[2968] | 320 | !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS |
---|
[3411] | 321 | CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:), & |
---|
| 322 | v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:))) |
---|
[2968] | 323 | |
---|
| 324 | !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV. |
---|
[3411] | 325 | i0=nbp_lev-locate(pintou(:),ptropou(i))+1 |
---|
| 326 | ll=check_ozone(v3(i, 1:i0 ,1),lon_in(i),lat_in(i),1 ,'troposphere', & |
---|
[2968] | 327 | 5.E-9,3.0E-7) |
---|
| 328 | ! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1) |
---|
[3411] | 329 | ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', & |
---|
[2968] | 330 | 5.E-8,1.5E-5) |
---|
| 331 | ! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1) |
---|
| 332 | |
---|
[2788] | 333 | END DO |
---|
| 334 | END IF |
---|
| 335 | |
---|
| 336 | CONTAINS |
---|
| 337 | |
---|
| 338 | |
---|
| 339 | !------------------------------------------------------------------------------- |
---|
| 340 | ! |
---|
| 341 | SUBROUTINE update_fields() |
---|
| 342 | ! |
---|
| 343 | !------------------------------------------------------------------------------- |
---|
| 344 | IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION |
---|
[2968] | 345 | CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.') |
---|
[2820] | 346 | irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number |
---|
| 347 | !--- MIN -> Security in the unlikely case of roundup errors. |
---|
[2788] | 348 | CALL get_3Dfields(v1) !--- Read ozone field(s) |
---|
[2820] | 349 | IF(lAdjTro) THEN !--- Additional files for fields strain |
---|
[3411] | 350 | IF(lPrSurf) CALL get_2Dfield(ps1,"ps") |
---|
| 351 | IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure") |
---|
| 352 | IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause") |
---|
[2788] | 353 | END IF |
---|
| 354 | ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS |
---|
[2981] | 355 | IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times |
---|
[2820] | 356 | IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN |
---|
[2981] | 357 | CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub) |
---|
[2788] | 358 | IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS |
---|
[2981] | 359 | WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '& |
---|
| 360 | //'reading record ',irec,' (time ',time_in(irec),')' |
---|
| 361 | CALL msg(.TRUE.,str,sub) |
---|
[2788] | 362 | CALL get_3Dfields(v1m) !--- Read ozone field(s) |
---|
[2820] | 363 | IF(lAdjTro) THEN !--- Additional files for fields strain |
---|
[3411] | 364 | IF(lPrSurf) CALL get_2Dfield(psm,"ps") |
---|
| 365 | IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure") |
---|
| 366 | IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause") |
---|
[2788] | 367 | END IF |
---|
| 368 | ELSE !=== SHIFT FIELDS |
---|
[2981] | 369 | irec=irec+1 |
---|
| 370 | WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'& |
---|
| 371 | //' current next one (',irec,', time ',time_in(irec),')' |
---|
| 372 | CALL msg(.TRUE.,str,sub) |
---|
[2788] | 373 | v1m=v1p !--- Ozone fields |
---|
[2820] | 374 | IF(lAdjTro) THEN !--- Additional files for fields strain |
---|
[3411] | 375 | IF(lPrSurf) psm=psp !--- Surface pressure |
---|
| 376 | IF(lPrTrop) ptm=ptp !--- Tropopause pressure |
---|
| 377 | IF(lO3Trop) otm=otp !--- Tropopause ozone |
---|
[2788] | 378 | END IF |
---|
| 379 | END IF |
---|
| 380 | irec=irec+1 |
---|
[2981] | 381 | WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'& |
---|
| 382 | ,irec,' (time ',time_in(irec),')' |
---|
| 383 | CALL msg(.TRUE.,str,sub) |
---|
[2788] | 384 | CALL get_3Dfields(v1p) !--- Read ozone field(s) |
---|
[2820] | 385 | IF(lAdjTro) THEN !--- Additional files for fields strain |
---|
[3411] | 386 | IF(lPrSurf) CALL get_2Dfield(psp,"ps") |
---|
| 387 | IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure") |
---|
| 388 | IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause") |
---|
[2788] | 389 | END IF |
---|
[2981] | 390 | irec=irec-1 |
---|
[2788] | 391 | END IF |
---|
| 392 | |
---|
| 393 | END SUBROUTINE update_fields |
---|
| 394 | ! |
---|
| 395 | !------------------------------------------------------------------------------- |
---|
| 396 | |
---|
| 397 | |
---|
| 398 | !------------------------------------------------------------------------------- |
---|
| 399 | ! |
---|
| 400 | SUBROUTINE get_2Dfield(v,var) |
---|
| 401 | ! |
---|
| 402 | !------------------------------------------------------------------------------- |
---|
| 403 | ! Purpose: Shortcut to get the "irec"th record of the surface field named "var" |
---|
| 404 | ! from the input file. |
---|
| 405 | ! Remark: In case the field is zonal, it is duplicated along first dimension. |
---|
| 406 | !------------------------------------------------------------------------------- |
---|
| 407 | ! Arguments: |
---|
| 408 | REAL, INTENT(INOUT) :: v(:,:) |
---|
| 409 | CHARACTER(LEN=*), INTENT(IN) :: var |
---|
| 410 | !------------------------------------------------------------------------------- |
---|
| 411 | CALL NF95_INQ_VARID(fID, TRIM(var), vID) |
---|
[2968] | 412 | CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim) |
---|
| 413 | IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[ 1,irec]) |
---|
| 414 | IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec]) |
---|
[2788] | 415 | CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID) |
---|
| 416 | |
---|
| 417 | !--- Flip latitudes: ascending in input file, descending in "rlatu". |
---|
[2968] | 418 | IF(n_dim==3) THEN |
---|
[3411] | 419 | v(1,:) = v(1,nbp_lat:1:-1) |
---|
| 420 | v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1) !--- Duplication |
---|
[2788] | 421 | ELSE |
---|
[3411] | 422 | v(:,:) = v(:,nbp_lat:1:-1) |
---|
[2788] | 423 | END IF |
---|
| 424 | |
---|
| 425 | END SUBROUTINE get_2Dfield |
---|
| 426 | ! |
---|
| 427 | !------------------------------------------------------------------------------- |
---|
| 428 | |
---|
| 429 | |
---|
| 430 | !------------------------------------------------------------------------------- |
---|
| 431 | ! |
---|
| 432 | SUBROUTINE get_3Dfields(v) |
---|
| 433 | ! |
---|
| 434 | !------------------------------------------------------------------------------- |
---|
| 435 | ! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam" |
---|
| 436 | ! from the input file. Fields are stacked on fourth dimension. |
---|
| 437 | ! Remark: In case the field is zonal, it is duplicated along first dimension. |
---|
| 438 | !------------------------------------------------------------------------------- |
---|
| 439 | ! Arguments: |
---|
| 440 | REAL, INTENT(INOUT) :: v(:,:,:,:) |
---|
| 441 | !------------------------------------------------------------------------------- |
---|
| 442 | DO i=1,SIZE(nam) |
---|
| 443 | CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID) |
---|
[2968] | 444 | CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim) |
---|
| 445 | IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[ 1,1,irec]) |
---|
| 446 | IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec]) |
---|
[2788] | 447 | CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID) |
---|
| 448 | END DO |
---|
| 449 | |
---|
| 450 | !--- Flip latitudes: ascending in input file, descending in "rlatu". |
---|
[2968] | 451 | IF(n_dim==3) THEN |
---|
[3411] | 452 | v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:) |
---|
| 453 | v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1) !--- Duplication |
---|
[2788] | 454 | ELSE |
---|
[3411] | 455 | v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:) |
---|
[2788] | 456 | END IF |
---|
| 457 | |
---|
| 458 | END SUBROUTINE get_3Dfields |
---|
| 459 | ! |
---|
| 460 | !------------------------------------------------------------------------------- |
---|
| 461 | |
---|
| 462 | |
---|
[2968] | 463 | |
---|
[2788] | 464 | !------------------------------------------------------------------------------- |
---|
| 465 | ! |
---|
[3411] | 466 | FUNCTION get_SigTrop(ih,it) |
---|
[2788] | 467 | ! |
---|
| 468 | !------------------------------------------------------------------------------- |
---|
| 469 | ! Arguments: |
---|
| 470 | INTEGER, INTENT(IN) :: ih |
---|
| 471 | INTEGER, INTENT(OUT) :: it |
---|
[3411] | 472 | REAL :: get_Sigtrop |
---|
[2788] | 473 | !------------------------------------------------------------------------------- |
---|
[3411] | 474 | !--- Pressure at tropopause is read in the forcing file |
---|
| 475 | IF(lPrTrop) THEN !--- PrTrop KNOWN FROM FILE |
---|
| 476 | get_SigTrop=pt2(ih)/ps2(ih); RETURN |
---|
| 477 | END IF |
---|
| 478 | !--- Chemical tropopause definition is used using a particular threshold |
---|
| 479 | IF(lO3Trop) THEN !--- o3trop KNOWN FROM FILE |
---|
| 480 | get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih)) |
---|
| 481 | ELSE IF(lo3tp) THEN !--- o3trop PARAMETRIZATION |
---|
| 482 | get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in) |
---|
| 483 | ELSE !--- o3trop CONSTANT |
---|
| 484 | get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0) |
---|
| 485 | END IF |
---|
| 486 | get_SigTrop=get_SigTrop/ps2(ih) |
---|
[2968] | 487 | |
---|
| 488 | END FUNCTION get_SigTrop |
---|
| 489 | ! |
---|
| 490 | !------------------------------------------------------------------------------- |
---|
| 491 | |
---|
| 492 | |
---|
| 493 | !------------------------------------------------------------------------------- |
---|
| 494 | ! |
---|
[3411] | 495 | FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop) |
---|
[2968] | 496 | ! |
---|
| 497 | !------------------------------------------------------------------------------- |
---|
| 498 | ! Purpose: Determine the tropopause using chemical criterium, ie the pressure |
---|
| 499 | ! at which the ozone concentration reaches a certain level. |
---|
| 500 | !------------------------------------------------------------------------------- |
---|
| 501 | ! Remarks: |
---|
| 502 | ! * Input field is upside down (increasing pressure // increasing vertical idx) |
---|
| 503 | ! The sweep is done from top to ground, starting at the 50hPa layer (idx it0), |
---|
| 504 | ! where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached: |
---|
| 505 | ! the O3 profile in this zone is decreasing with pressure. |
---|
| 506 | ! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or |
---|
| 507 | ! determined using an empirical law roughly derived from ... & al. |
---|
| 508 | !------------------------------------------------------------------------------- |
---|
| 509 | ! Arguments: |
---|
[3411] | 510 | REAL :: chem_tropopause !--- Pressure at tropopause |
---|
[2968] | 511 | INTEGER, INTENT(IN) :: ih !--- Horizontal index |
---|
| 512 | INTEGER, INTENT(OUT) :: it !--- Index of tropopause layer |
---|
| 513 | INTEGER, INTENT(IN) :: it0 !--- Idx: higher than tropopause |
---|
[3411] | 514 | REAL, INTENT(IN) :: pint(:) !--- Cells-interf Pr, increasing |
---|
[2968] | 515 | REAL, INTENT(IN) :: o3(:,:) !--- Ozone field (pptV) |
---|
[3411] | 516 | REAL, OPTIONAL, INTENT(IN) :: pcen(:) !--- Cells-center Pr, increasing |
---|
[2968] | 517 | REAL, OPTIONAL, INTENT(IN) :: o3trop !--- Ozone at tropopause |
---|
| 518 | !------------------------------------------------------------------------------- |
---|
[2788] | 519 | ! Local variables: |
---|
[3411] | 520 | REAL :: o3t !--- Ozone concent. at tropopause |
---|
| 521 | REAL :: al !--- Interpolation coefficient |
---|
| 522 | REAL :: coef !--- Coeff of latitude modulation |
---|
[2968] | 523 | REAL, PARAMETER :: co3(3)=[91.,28.,20.] !--- Coeff for o3 at tropopause |
---|
[2788] | 524 | !------------------------------------------------------------------------------- |
---|
[2968] | 525 | !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN |
---|
| 526 | IF(PRESENT(o3trop)) THEN !=== THRESHOLD FROM ARGUMENTS |
---|
| 527 | o3t=o3trop |
---|
| 528 | ELSE !=== EMPIRICAL LAW |
---|
| 529 | coef = TANH(lat_in(ih)/co3(3)) !--- Latitude modulation |
---|
| 530 | coef = SIN (2.*RPI*(y_frac-1./6.)) * coef !--- Seasonal modulation |
---|
[2788] | 531 | o3t = 1.E-9 * (co3(1)+co3(2)*coef) !--- Value from parametrization |
---|
| 532 | END IF |
---|
| 533 | |
---|
[2968] | 534 | !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1) |
---|
| 535 | it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO |
---|
| 536 | al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1)) |
---|
[3411] | 537 | IF(PRESENT(pcen)) THEN |
---|
| 538 | chem_tropopause = pcen(it)**(1.-al) * pcen(it+1)**al |
---|
| 539 | ELSE |
---|
| 540 | chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) ) |
---|
| 541 | END IF |
---|
| 542 | it = locate(pint(:), chem_tropopause) !--- pint(it)<ptrop<pint(it+1) |
---|
[2968] | 543 | |
---|
[3411] | 544 | END FUNCTION chem_tropopause |
---|
[2788] | 545 | ! |
---|
| 546 | !------------------------------------------------------------------------------- |
---|
| 547 | |
---|
| 548 | |
---|
[2968] | 549 | !------------------------------------------------------------------------------- |
---|
| 550 | ! |
---|
[3411] | 551 | FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) |
---|
[2968] | 552 | ! |
---|
| 553 | !------------------------------------------------------------------------------- |
---|
| 554 | IMPLICIT NONE |
---|
| 555 | !------------------------------------------------------------------------------- |
---|
| 556 | ! Arguments: |
---|
[3411] | 557 | LOGICAL :: check_ozone !--- .T. => some wrong values |
---|
[2968] | 558 | REAL, INTENT(IN) :: o3col(:), lon, lat |
---|
| 559 | INTEGER, INTENT(IN) :: ilev0 |
---|
| 560 | CHARACTER(LEN=*), INTENT(IN) :: layer |
---|
| 561 | REAL, OPTIONAL, INTENT(IN) :: vmin |
---|
| 562 | REAL, OPTIONAL, INTENT(IN) :: vmax |
---|
| 563 | !------------------------------------------------------------------------------- |
---|
| 564 | ! Local variables: |
---|
| 565 | INTEGER :: k |
---|
| 566 | LOGICAL :: lmin, lmax |
---|
| 567 | REAL :: fac |
---|
| 568 | CHARACTER(LEN=6) :: unt |
---|
| 569 | !------------------------------------------------------------------------------- |
---|
| 570 | !--- NOTHING TO DO |
---|
| 571 | lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0 |
---|
| 572 | lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0 |
---|
[3411] | 573 | check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN |
---|
[2968] | 574 | |
---|
| 575 | !--- SOME TOO LOW VALUES FOUND |
---|
| 576 | IF(lmin) THEN |
---|
| 577 | CALL unitc(vmin,unt,fac) |
---|
| 578 | DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE |
---|
| 579 | WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' & |
---|
| 580 | //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', & |
---|
| 581 | fac*vmin,unt//' in '//TRIM(layer) |
---|
| 582 | END DO |
---|
| 583 | END IF |
---|
| 584 | |
---|
| 585 | !--- SOME TOO HIGH VALUES FOUND |
---|
| 586 | IF(lmax) THEN |
---|
| 587 | CALL unitc(vmax,unt,fac) |
---|
| 588 | DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE |
---|
| 589 | WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' & |
---|
| 590 | //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', & |
---|
| 591 | fac*vmax,unt//' in '//TRIM(layer) |
---|
| 592 | END DO |
---|
| 593 | END IF |
---|
| 594 | |
---|
| 595 | END FUNCTION check_ozone |
---|
| 596 | ! |
---|
| 597 | !------------------------------------------------------------------------------- |
---|
| 598 | |
---|
| 599 | |
---|
| 600 | !------------------------------------------------------------------------------- |
---|
| 601 | ! |
---|
| 602 | SUBROUTINE unitc(val,unt,fac) |
---|
| 603 | ! |
---|
| 604 | !------------------------------------------------------------------------------- |
---|
| 605 | IMPLICIT NONE |
---|
| 606 | !------------------------------------------------------------------------------- |
---|
| 607 | ! Arguments: |
---|
| 608 | REAL, INTENT(IN) :: val |
---|
| 609 | CHARACTER(LEN=6), INTENT(OUT) :: unt |
---|
| 610 | REAL, INTENT(OUT) :: fac |
---|
| 611 | !------------------------------------------------------------------------------- |
---|
| 612 | ! Local variables: |
---|
| 613 | INTEGER :: ndgt |
---|
| 614 | !------------------------------------------------------------------------------- |
---|
| 615 | ndgt=3*FLOOR(LOG10(val)/3.) |
---|
| 616 | SELECT CASE(ndgt) |
---|
| 617 | CASE(-6); unt=' ppmV '; fac=1.E6 |
---|
| 618 | CASE(-9); unt=' ppbV '; fac=1.E9 |
---|
| 619 | CASE DEFAULT; unt=' vmr '; fac=1.0 |
---|
| 620 | END SELECT |
---|
| 621 | |
---|
| 622 | END SUBROUTINE unitc |
---|
| 623 | ! |
---|
| 624 | !------------------------------------------------------------------------------- |
---|
| 625 | |
---|
| 626 | |
---|
| 627 | !------------------------------------------------------------------------------- |
---|
| 628 | ! |
---|
| 629 | SUBROUTINE msg(ll,str,sub) |
---|
| 630 | ! |
---|
| 631 | !------------------------------------------------------------------------------- |
---|
| 632 | USE print_control_mod, ONLY: lunout |
---|
| 633 | IMPLICIT NONE |
---|
| 634 | !------------------------------------------------------------------------------- |
---|
| 635 | ! Arguments: |
---|
| 636 | LOGICAL, INTENT(IN) :: ll |
---|
| 637 | CHARACTER(LEN=*), INTENT(IN) :: str |
---|
| 638 | CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub |
---|
| 639 | !------------------------------------------------------------------------------- |
---|
| 640 | IF(.NOT.ll) RETURN |
---|
| 641 | IF(PRESENT(sub)) THEN |
---|
| 642 | WRITE(lunout,*)TRIM(sub)//': '//TRIM(str) |
---|
| 643 | ELSE |
---|
| 644 | WRITE(lunout,*)TRIM(str) |
---|
| 645 | END IF |
---|
| 646 | |
---|
| 647 | END SUBROUTINE msg |
---|
| 648 | ! |
---|
| 649 | !------------------------------------------------------------------------------- |
---|
| 650 | |
---|
[2788] | 651 | END SUBROUTINE regr_pr_time_av |
---|
| 652 | ! |
---|
| 653 | !------------------------------------------------------------------------------- |
---|
| 654 | |
---|
| 655 | END MODULE regr_pr_time_av_m |
---|