! $Id$ MODULE regr_pr_time_av_m USE write_field_phy USE mod_phys_lmdz_transfert_para, ONLY: bcast USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank USE print_control_mod, ONLY: prt_level IMPLICIT NONE !------------------------------------------------------------------------------- ! Purpose: Regrid pressure with an averaging method. Operations done: ! * on the root process: read a NetCDF 3D or 4D field at a given day. ! * pack the fields to the LMDZ "horizontal "physics" grid and scatter ! to all threads of all processes; ! * in all the threads of all the processes, regrid the fields in pressure ! to the LMDZ vertical grid. ! * the forcing fields are stretched if the following arguments are present: ! - "lat_in": input file latitudes. ! - "pcen_in": input file cells center pressure levels. ! - "ptrop_ou": target grid (LMDZ) tropopause pressure. ! so that the tropopause is shifted to the position of the one of LMDZ. ! Note that the ozone quantity conservation is not ensured. !------------------------------------------------------------------------------- ! Initial routine: regr_pr_av_m module (L. Guez). ! Present version: David Cugnet ; corresponding additions: ! - time interpolation ! - 3D compliant ! - vertical stretching of the field to allow tropopause matching between ! input field and current lmdz field. !------------------------------------------------------------------------------- ! Remarks: ! * 3D fields are zonal means, with dimensions (latitude, pressure, julian day) ! * 4D fields have the dimensions: (longitude, latitude, pressure, julian day) ! * All the fields are already on the horizontal grid (rlatu) or (rlatu,rlonv), ! except that the latitudes are in ascending order in the input file. ! * We assume that all the input fields have the same coordinates. ! * The target vertical LMDZ grid is the grid of layer boundaries. ! * Regridding in pressure is conservative, second order. ! * All the fields are regridded as a single multi-dimensional array, so it ! saves CPU time to call this procedure once for several NetCDF variables ! rather than several times, each time for a single NetCDF variable. ! * The input file pressure at tropopause can be (in decreasing priority): ! 1) read from the file if "tropopause_air_pressure" field is available. ! 2) computed using "tro3" and "tro3_at_tropopause' (if available). ! 3) computed using "tro3" and a fixed threshold otherwise, determined using ! an empirical three parameters law: ! o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3) ! => co1 and co2 are in ppbV, and co3 in degrees. ! => co3 allow a smooth transition between north and south hemispheres. ! * If available, the field "ps" (input file pressure at surface) is used. ! If not, the current LMDZ ground pressure is taken instead. ! * Fields with suffix "m"/"p" are at the closest records earlier/later than ! the mid-julian day "julien", on the global "dynamics" horizontal grid. ! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure ! interval "pint_in(k:k+1)]" and variable "nam(l)". ! * In the 2D file case, the values are the same for all longitudes. ! * The tropopause correction works like this: the input fields (file) are ! interpolated on output (LMDZ) pressure field, which is streched using a power ! law in a limited zone made of 3 layers: ! 1) between the two tropopauses (file and LMDZ) ! 2) in an upper and a lower transitions layers ! The stretching function has the following form: ! Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou) ! This value shifts the file tropopause to the height of the one of LMDZ. ! The stretching is fully applied in the central zone only, and only partially ! in the transitions zones, thick enough to guarantee a growing stretched ! pressure field. The ponderation function for alpha to modulate the stretching ! is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0) ! in the transition layers (from 1 to 0 ; in fact: sections of 1/log function), ! making this localisation function quasi-trapezoidal. ! ! * The following fields are on the global "dynamics" grid, as read from files: REAL, SAVE, ALLOCATABLE :: v1m(:,:,:,:) !--- Previous ozone fields REAL, SAVE, ALLOCATABLE :: v1p(:,:,:,:) !--- Next ozone fields REAL, SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure REAL, SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure REAL, SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio INTEGER, SAVE :: ntim_in !--- Records nb in input file INTEGER, SAVE :: itrp0 !--- idx above chem tropop. INTEGER, SAVE :: irec !--- Current time index ! * for daily input files: current julian day number ! * for monthly input files: julien is in [time_in(irec),time_in(irec+1)] LOGICAL, SAVE :: linterp !--- Interpolation in time LOGICAL, SAVE :: lPrSurf !--- Surface pressure flag LOGICAL, SAVE :: lPrTrop !--- Tropopause pressure flag LOGICAL, SAVE :: lO3Trop !--- Tropopause ozone flag LOGICAL, SAVE :: lfirst=.TRUE. !--- First call flag !$OMP THREADPRIVATE(lfirst) REAL, PARAMETER :: pTropUp=9.E+3 !--- Value < tropopause pressure (Pa) REAL, PARAMETER :: gamm = 0.4 !--- Relative thickness of transitions REAL, PARAMETER :: rho = 1.3 !--- Max tropopauses sigma ratio REAL, PARAMETER :: o3t0 = 1.E-7 !--- Nominal O3 vmr at tropopause LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause LOGICAL, PARAMETER :: debug=.FALSE. !--- Force writefield_phy multiple outputs REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum REAL, PARAMETER :: ChemPTrMax=4.E+4 ! chemical tropopause pressure (Pa). REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum REAL, PARAMETER :: DynPTrMax =4.E+4 ! dynamical tropopause pressure (Pa). CONTAINS !------------------------------------------------------------------------------- ! SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, & time_in, lon_in, lat_in, pcen_in, ptrop_ou) ! !------------------------------------------------------------------------------- USE dimphy, ONLY: klon USE netcdf95, ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, & NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION USE netcdf, ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR USE assert_m, ONLY: assert USE assert_eq_m, ONLY: assert_eq USE comvert_mod, ONLY: scaleheight USE interpolation, ONLY: locate USE regr_conserv_m, ONLY: regr_conserv USE regr_lint_m, ONLY: regr_lint USE slopes_m, ONLY: slopes USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter USE phys_cal_mod, ONLY: calend, year_len, days_elapsed, jH_cur !------------------------------------------------------------------------------- ! Arguments: INTEGER, INTENT(IN) :: fID !--- NetCDF file ID CHARACTER(LEN=13), INTENT(IN) :: nam(:) !--- NetCDF variables names REAL, INTENT(IN) :: julien !--- Days since Jan 1st !--- File & LMDZ (resp. decreasing & increasing order) interfaces pressure, Pa REAL, INTENT(IN) :: pint_in(:) !--- in: file (nlev_in+1) REAL, INTENT(IN) :: pint_ou(:,:) !--- out: LMDZ (klon,nbp_lev+1) REAL, INTENT(OUT) :: v3(:,:,:) !--- Regridded fld (klon,llm,n_var) REAL, INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days ! since Jan 1 of current year REAL, INTENT(IN), OPTIONAL :: lon_in(:) !--- File longitudes vector (klon) REAL, INTENT(IN), OPTIONAL :: lat_in(:) !--- File latitudes vector (klon) REAL, INTENT(IN), OPTIONAL :: pcen_in(:) !--- File mid-layers pres (nlev_in) REAL, INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- LMDZ tropopause pres (klon) !------------------------------------------------------------------------------- ! Local variables: include "clesphys.h" include "YOMCST.h" CHARACTER(LEN=80) :: sub CHARACTER(LEN=320) :: str INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp INTEGER :: i, nlev_in, n_dim, itop, ito0, i0 LOGICAL :: lAdjTro !--- Need to adjust tropopause REAL :: y_frac !--- Elapsed year fraction REAL :: alpha, beta, al !--- For stretching/interpolation REAL :: SigT_in, SigT_ou !--- Input and output tropopauses REAL :: Sig_bot, Sig_top !--- Fully strained layer bounds REAL :: Sig_bo0, Sig_to0 !--- Total strained layers bounds REAL :: Sig_in(SIZE(pint_in)) !--- Input field sigma levels REAL :: Sig_ou (nbp_lev+1) !--- Output LMDZ sigma levels REAL :: phi (nbp_lev+1) !--- Stretching exponent anomaly REAL :: pstr_ou(nbp_lev+1) !--- Stretched pressure levels REAL :: pintou (nbp_lev+1) !--- pint_ou in reversed order REAL :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam)) REAL :: v2(klon, SIZE(pint_in)-1,SIZE(nam)) ! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid. ! v2: Field scattered to the partial "physics" horizontal grid. ! In the 2D file case, the values are the same for all longitudes. ! "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)". ! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)" REAL, DIMENSION(nbp_lon, nbp_lat) :: ps1, pt1, ot1 REAL, DIMENSION(klon) :: ps2, pt2, ot2, ptropou LOGICAL :: ll !--- For debug REAL, DIMENSION(klon) :: ptrop_in, ptrop_ef!, dzStrain, dzStrain0 ! REAL, DIMENSION(klon,nbp_lev+1) :: pstrn_ou, phii !------------------------------------------------------------------------------- sub="regr_pr_time_av" nlev_in=SIZE(pint_in)-1 CALL assert(SIZE(v3,1)==klon, TRIM(sub)//" v3 klon") CALL assert(SIZE(v3,2)==nbp_lev, TRIM(sub)//" v3 nbp_lev") n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var") CALL assert(SIZE(pint_ou,1)==klon ,TRIM(sub)//" pint_ou klon") CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1") IF(PRESENT(lon_in)) CALL assert(SIZE(lon_in )==klon,TRIM(sub)//" lon_in klon") IF(PRESENT(lat_in)) CALL assert(SIZE(lat_in )==klon,TRIM(sub)//" lat_in klon") IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon") IF(PRESENT(pcen_in)) CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in") lAdjTro=PRESENT(ptrop_ou) IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) & CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1) !$OMP MASTER IF(is_mpi_root) THEN !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE IF(lfirst) THEN lPrSurf=NF90_INQ_VARID(fID,"ps" ,vID)==NF90_NOERR lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause" ,vID)==NF90_NOERR CALL NF95_INQ_DIMID(fID,"time",vID) CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in) linterp=PRESENT(time_in).AND.ntim_in==14 IF(linterp) THEN ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var)) ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var)) ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat)) ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat)) IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat)) END IF !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa) IF(lAdjTro) itrp0=locate(pcen_in,pTropUp) CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub) CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub) CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub) IF(lPrTrop) THEN CALL msg(lAdjTro,' PRESSURE AT TROPOPAUSE from file.') ELSE IF(lO3Trop) THEN CALL msg(lAdjTro,' O3 CONCENTRATION AT TROPOPAUSE from file.') ELSE IF(lo3tp) THEN CALL msg(lAdjTro,' PARAMETRIZED O3 concentration at tropopause.') ELSE CALL msg(lAdjTro,' CONSTANT O3 concentration at tropopause.') END IF END IF !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES) CALL update_fields() !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES IF(linterp) THEN WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '& &//'fields at times ',time_in(irec),' and ', time_in(irec+1) CALL msg(.TRUE.,str,sub) al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec)) v1=al*v1m+(1.-al)*v1p IF(lPrSurf) ps1=al*psm+(1.-al)*psp IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp IF(lO3Trop) ot1=al*otm+(1.-al)*otp END IF END IF !$OMP END MASTER IF(lfirst) THEN lfirst=.FALSE.; CALL bcast(lfirst) IF(lAdjTro) CALL bcast(itrp0) CALL bcast(lPrTrop); CALL bcast(lPrSurf) CALL bcast(lO3Trop); CALL bcast(linterp) END IF CALL scatter2d(v1,v2) !--- No "ps" in input file => assumed to be equal to current LMDZ ground press IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF IF(lPrTrop) CALL scatter2d(pt1,pt2) IF(lO3Trop) CALL scatter2d(ot1,ot2) !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing IF(.NOT.lAdjTro) THEN DO i=1,klon pintou = pint_ou(i,nbp_lev+1:1:-1) CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:), & v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:))) END DO ELSE y_frac=(REAL(days_elapsed)+jH_cur)/year_len !--- OUTPUT SIGMA LEVELS DO i=1,klon !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES pintou=pint_ou(i,nbp_lev+1:1:-1) !--- increasing values Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)] !--- increasing values Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered ! to keep tropopause pressure realistic ; high values are usually due to ! ozone hole fooling the crude chemical tropopause detection algorithm. SigT_in = get_SigTrop(i,itrp) SigT_in=MIN(SigT_in,ChemPTrMax/ps2(i)) !--- too low value filtered SigT_in=MAX(SigT_in,ChemPTrMin/ps2(i)) !--- too high value filtered !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the ! dynamical tropopause (especially in filaments) are conterbalanced with a ! filter ensuring it stays within a certain distance around input (file) ! tropopause, hence avoiding avoid a too thick stretched region ; a final ! extra-safety filter keeps the tropopause pressure value realistic. SigT_ou = ptrop_ou(i)/ps2(i) IF(SigT_ouSigT_in*rho) SigT_ou=SigT_in*rho !--- too high value w/r input SigT_ou=MIN(SigT_ou,DynPTrMax/ps2(i)) !--- too low value filtered SigT_ou=MAX(SigT_ou,DynPTrMin/ps2(i)) !--- too high value filtered ptropou(i)=SigT_ou*ps2(i) !--- FULLY STRETCHED LAYER BOUNDS (alpha power law fully applied) alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou) Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot) Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top) !--- PARTIALLY STRETCHED LAYER BOUNDS (fraction of alpha applied) ! The used formulae ensure the first derivative stays positive. beta = LOG(Sig_top)/LOG(Sig_bot) Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta) Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top ** beta !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS !--- gamma1 in [Sig_to0,Sig_top] 1->0 in [Sig_bot,Sig_bo0] phi(:)=0. phi(itop+1:ibot) = 1. phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))& *LOG(Sig_top)/LOG(Sig_top/Sig_to0) phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))& *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0) !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER) pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:)) !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:), & v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:))) !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV. i0=nbp_lev-locate(pintou(:),ptropou(i))+1 ll=check_ozone(v3(i, 1:i0 ,1),lon_in(i),lat_in(i),1 ,'troposphere', & 5.E-9,3.0E-7) ! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1) ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', & 5.E-8,1.5E-5) ! IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1) IF(debug) THEN ! dzStrain0(i)=7.*LOG(Sig_bo0/Sig_to0) ! dzStrain (i)=7.*LOG(Sig_bot/Sig_top) ptrop_in(i) = SigT_in*ps2(i) ! Pstrn_ou(i,:)= pstr_ou ! phii(i,:)=phi(:) ptrop_ef(i)=chem_tropopause(i, itrp, locate(pintou(:),PTropUp), & pintou(:), v3(:,nbp_lev+1:1:-1,1),o3trop=o3t0) END IF END DO IF(debug) THEN ! CALL writefield_phy('PreSt_ou' ,Pstrn_ou,nbp_lev+1) !--- Strained pressure ! CALL writefield_phy('dzStrain' ,dzStrain ,1) !--- Fully & total strained ! CALL writefield_phy('dzStrain0',dzStrain0,1) !--- regions thickness ! CALL writefield_phy('phi',phii,nbp_lev+1) !--- Localisation function !--- Tropopauses pressures: CALL writefield_phy('PreTr_in',ptrop_in,1) !--- Input and effective CALL writefield_phy('PreTr_ef',ptrop_ef,1) ! chemical tropopauses END IF END IF IF(debug) THEN IF(lAdjTro) & CALL writefield_phy('PreTr_ou',ptropou,1) !--- LMDz dyn tropopause CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field CALL writefield_phy('Ozone_ou',v3(:,:,1),nbp_lev)!--- Output ozone field CALL writefield_phy('Pint_ou' ,pint_ou,1+nbp_lev)!--- LMDZ unstreched Press END IF CONTAINS !------------------------------------------------------------------------------- ! SUBROUTINE update_fields() ! !------------------------------------------------------------------------------- IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.') irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number !--- MIN -> Security in the unlikely case of roundup errors. CALL get_3Dfields(v1) !--- Read ozone field(s) IF(lAdjTro) THEN !--- Additional files for fields strain IF(lPrSurf) CALL get_2Dfield(ps1,"ps") IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure") IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause") END IF ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times IF(.NOT.lfirst.AND.julieno3t>o3(ih,it+1) it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1)) IF(PRESENT(pcen)) THEN chem_tropopause = pcen(it)**(1.-al) * pcen(it+1)**al ELSE chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) ) END IF it = locate(pint(:), chem_tropopause) !--- pint(it) some wrong values REAL, INTENT(IN) :: o3col(:), lon, lat INTEGER, INTENT(IN) :: ilev0 CHARACTER(LEN=*), INTENT(IN) :: layer REAL, OPTIONAL, INTENT(IN) :: vmin REAL, OPTIONAL, INTENT(IN) :: vmax !------------------------------------------------------------------------------- ! Local variables: INTEGER :: k LOGICAL :: lmin, lmax REAL :: fac CHARACTER(LEN=6) :: unt !------------------------------------------------------------------------------- !--- NOTHING TO DO lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3colvmax)/=0 check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN IF(prt_level<1) RETURN !--- SOME TOO LOW VALUES FOUND IF(lmin) THEN CALL unitc(vmin,unt,fac) DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in ' & //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', & fac*vmin,unt//' in '//TRIM(layer) END DO END IF !--- SOME TOO HIGH VALUES FOUND IF(lmax) THEN CALL unitc(vmax,unt,fac) DO k=1,SIZE(o3col); IF(o3col(k) ', & fac*vmax,unt//' in '//TRIM(layer) END DO END IF END FUNCTION check_ozone ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE unitc(val,unt,fac) ! !------------------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: REAL, INTENT(IN) :: val CHARACTER(LEN=6), INTENT(OUT) :: unt REAL, INTENT(OUT) :: fac !------------------------------------------------------------------------------- ! Local variables: INTEGER :: ndgt !------------------------------------------------------------------------------- ndgt=3*FLOOR(LOG10(val)/3.) SELECT CASE(ndgt) CASE(-6); unt=' ppmV '; fac=1.E6 CASE(-9); unt=' ppbV '; fac=1.E9 CASE DEFAULT; unt=' vmr '; fac=1.0 END SELECT END SUBROUTINE unitc ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE msg(ll,str,sub) ! !------------------------------------------------------------------------------- USE print_control_mod, ONLY: lunout IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: LOGICAL, INTENT(IN) :: ll CHARACTER(LEN=*), INTENT(IN) :: str CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub !------------------------------------------------------------------------------- IF(.NOT.ll) RETURN IF(PRESENT(sub)) THEN WRITE(lunout,*)TRIM(sub)//': '//TRIM(str) ELSE WRITE(lunout,*)TRIM(str) END IF END SUBROUTINE msg ! !------------------------------------------------------------------------------- END SUBROUTINE regr_pr_time_av ! !------------------------------------------------------------------------------- END MODULE regr_pr_time_av_m