MODULE diaobs !!====================================================================== !! *** MODULE diaobs *** !! Observation diagnostics: Computation of the misfit between data and !! their model equivalent !!====================================================================== !!---------------------------------------------------------------------- !! dia_obs_init : Reading and prepare observations !! dia_obs : Compute model equivalent to observations !! dia_obs_wri : Write observational diagnostics !! ini_date : Compute the initial date YYYYMMDD.HHMMSS !! fin_date : Compute the final date YYYYMMDD.HHMMSS !!---------------------------------------------------------------------- !! * Modules used USE wrk_nemo ! Memory Allocation USE par_kind ! Precision variables USE in_out_manager ! I/O manager USE par_oce USE dom_oce ! Ocean space and time domain variables USE obs_read_prof ! Reading and allocation of profile obs USE obs_read_surf ! Reading and allocation of surface obs USE obs_readmdt ! Reading and allocation of MDT for SLA. USE obs_prep ! Preparation of obs. (grid search etc). USE obs_oper ! Observation operators USE obs_write ! Writing of observation related diagnostics USE obs_grid ! Grid searching USE obs_read_altbias ! Bias treatment for altimeter USE obs_sstbias ! Bias correction routine for SST USE obs_profiles_def ! Profile data definitions USE obs_surf_def ! Surface data definitions USE obs_types ! Definitions for observation types USE mpp_map ! MPP mapping USE lib_mpp ! For ctl_warn/stop IMPLICIT NONE !! * Routine accessibility PRIVATE PUBLIC dia_obs_init, & ! Initialize and read observations & dia_obs, & ! Compute model equivalent to observations & dia_obs_wri, & ! Write model equivalent to observations & dia_obs_dealloc ! Deallocate dia_obs data !! * Module variables LOGICAL, PUBLIC :: & & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. LOGICAL :: ln_diaobs !: Logical switch for the obs operator LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint INTEGER :: nn_1dint !: Vertical interpolation method INTEGER :: nn_2dint_default !: Default horizontal interpolation method INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (-1 = default) INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (-1 = default) INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (-1 = default) INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method (-1 = default) INTEGER, DIMENSION(imaxavtypes) :: & & nn_profdavtypes !: Profile data types representing a daily average INTEGER :: nproftypes !: Number of profile obs types INTEGER :: nsurftypes !: Number of surface obs types INTEGER, DIMENSION(:), ALLOCATABLE :: & & nvarsprof, & !: Number of profile variables & nvarssurf !: Number of surface variables INTEGER, DIMENSION(:), ALLOCATABLE :: & & nextrprof, & !: Number of profile extra variables & nextrsurf !: Number of surface extra variables INTEGER, DIMENSION(:), ALLOCATABLE :: & & n2dintsurf !: Interpolation option for surface variables REAL(wp), DIMENSION(:), ALLOCATABLE :: & & ravglamscl, & !: E/W diameter of averaging footprint for surface variables & ravgphiscl !: N/S diameter of averaging footprint for surface variables LOGICAL, DIMENSION(:), ALLOCATABLE :: & & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres & llnightav !: Logical for calculating night-time averages TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & & surfdata, & !: Initial surface data & surfdataqc !: Surface data after quality control TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & & profdata, & !: Initial profile data & profdataqc !: Profile data after quality control CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & & cobstypesprof, & !: Profile obs types & cobstypessurf !: Surface obs types !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- !! * Substitutions # include "domzgr_substitute.h90" CONTAINS SUBROUTINE dia_obs_init !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_init *** !! !! ** Purpose : Initialize and read observations !! !! ** Method : Read the namelist and call reading routines !! !! ** Action : Read the namelist and call reading routines !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (A. Weaver) Reformatted !! ! 06-10 (A. Weaver) Cleaning and add controls !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 14-08 (J.While) Incorporated SST bias correction !! ! 15-02 (M. Martin) Simplification of namelist and code !!---------------------------------------------------------------------- IMPLICIT NONE !! * Local declarations INTEGER, PARAMETER :: & & jpmaxnfiles = 1000 ! Maximum number of files for each obs type INTEGER, DIMENSION(:), ALLOCATABLE :: & & ifilesprof, & ! Number of profile files & ifilessurf ! Number of surface files INTEGER :: ios ! Local integer output status for namelist read INTEGER :: jtype ! Counter for obs types INTEGER :: jvar ! Counter for variables INTEGER :: jfile ! Counter for files INTEGER :: jnumsstbias ! Number of SST bias files to read and apply INTEGER :: n2dint_type ! Local version of nn_2dint* CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & & cn_profbfiles, & ! T/S profile input filenames & cn_sstfbfiles, & ! Sea surface temperature input filenames & cn_slafbfiles, & ! Sea level anomaly input filenames & cn_sicfbfiles, & ! Seaice concentration input filenames & cn_velfbfiles, & ! Velocity profile input filenames & cn_sssfbfiles, & ! Sea surface salinity input filenames & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames & cn_pno3fbfiles, & ! Profile nitrate input filenames & cn_psi4fbfiles, & ! Profile silicate input filenames & cn_ppo4fbfiles, & ! Profile phosphate input filenames & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames & cn_palkfbfiles, & ! Profile alkalinity input filenames & cn_pphfbfiles, & ! Profile pH input filenames & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames & cn_sstbiasfiles ! SST bias input filenames CHARACTER(LEN=128) :: & & cn_altbiasfile ! Altimeter bias input filename LOGICAL :: ln_t3d ! Logical switch for temperature profiles LOGICAL :: ln_s3d ! Logical switch for salinity profiles LOGICAL :: ln_sla ! Logical switch for sea level anomalies LOGICAL :: ln_sst ! Logical switch for sea surface temperature LOGICAL :: ln_sic ! Logical switch for sea ice concentration LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs LOGICAL :: ln_slchlnon ! Logical switch for surface non-diatom log10(chlorophyll) obs LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs LOGICAL :: ln_slphynon ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs LOGICAL :: ln_pph ! Logical switch for profile pH obs LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs LOGICAL :: ln_nea ! Logical switch to remove obs near land LOGICAL :: ln_altbias ! Logical switch for altimeter bias LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & & clproffiles, & ! Profile filenames & clsurffiles ! Surface filenames LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zglam ! Model longitudes for profile variables REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zgphi ! Model latitudes for profile variables REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & & zmask ! Model land/sea mask associated with variables NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & & ln_sst, ln_sic, ln_sss, ln_vel3d, & & ln_slchltot, ln_slchldia, ln_slchlnon, & & ln_slchldin, ln_slchlmic, ln_slchlnan, & & ln_slchlpic, ln_schltot, & & ln_slphytot, ln_slphydia, ln_slphynon, & & ln_sspm, ln_sfco2, ln_spco2, & & ln_plchltot, ln_pchltot, ln_pno3, & & ln_psi4, ln_ppo4, ln_pdic, & & ln_palk, ln_pph, ln_po2, & & ln_altbias, ln_sstbias, ln_nea, & & ln_grid_global, ln_grid_search_lookup, & & ln_ignmis, ln_s_at_t, ln_bound_reject, & & ln_sstnight, ln_default_fp_indegs, & & ln_sla_fp_indegs, ln_sst_fp_indegs, & & ln_sss_fp_indegs, ln_sic_fp_indegs, & & cn_profbfiles, cn_slafbfiles, & & cn_sstfbfiles, cn_sicfbfiles, & & cn_velfbfiles, cn_sssfbfiles, & & cn_slchltotfbfiles, cn_slchldiafbfiles, & & cn_slchlnonfbfiles, cn_slchldinfbfiles, & & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & & cn_slchlpicfbfiles, cn_schltotfbfiles, & & cn_slphytotfbfiles, cn_slphydiafbfiles, & & cn_slphynonfbfiles, cn_sspmfbfiles, & & cn_sfco2fbfiles, cn_spco2fbfiles, & & cn_plchltotfbfiles, cn_pchltotfbfiles, & & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & & cn_po2fbfiles, & & cn_sstbiasfiles, cn_altbiasfile, & & cn_gridsearchfile, rn_gridsearchres, & & rn_dobsini, rn_dobsend, & & rn_default_avglamscl, rn_default_avgphiscl, & & rn_sla_avglamscl, rn_sla_avgphiscl, & & rn_sst_avglamscl, rn_sst_avgphiscl, & & rn_sss_avglamscl, rn_sss_avgphiscl, & & rn_sic_avglamscl, rn_sic_avgphiscl, & & nn_1dint, nn_2dint_default, & & nn_2dint_sla, nn_2dint_sst, & & nn_2dint_sss, nn_2dint_sic, & & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & & nn_profdavtypes !----------------------------------------------------------------------- ! Read namelist parameters !----------------------------------------------------------------------- ! Some namelist arrays need initialising cn_profbfiles(:) = '' cn_slafbfiles(:) = '' cn_sstfbfiles(:) = '' cn_sicfbfiles(:) = '' cn_velfbfiles(:) = '' cn_sssfbfiles(:) = '' cn_slchltotfbfiles(:) = '' cn_slchldiafbfiles(:) = '' cn_slchlnonfbfiles(:) = '' cn_slchldinfbfiles(:) = '' cn_slchlmicfbfiles(:) = '' cn_slchlnanfbfiles(:) = '' cn_slchlpicfbfiles(:) = '' cn_schltotfbfiles(:) = '' cn_slphytotfbfiles(:) = '' cn_slphydiafbfiles(:) = '' cn_slphynonfbfiles(:) = '' cn_sspmfbfiles(:) = '' cn_sfco2fbfiles(:) = '' cn_spco2fbfiles(:) = '' cn_plchltotfbfiles(:) = '' cn_pchltotfbfiles(:) = '' cn_pno3fbfiles(:) = '' cn_psi4fbfiles(:) = '' cn_ppo4fbfiles(:) = '' cn_pdicfbfiles(:) = '' cn_palkfbfiles(:) = '' cn_pphfbfiles(:) = '' cn_po2fbfiles(:) = '' cn_sstbiasfiles(:) = '' nn_profdavtypes(:) = -1 CALL ini_date( rn_dobsini ) CALL fin_date( rn_dobsend ) ! Read namelist namobs : control observation diagnostics REWIND( numnam_ref ) ! Namelist namobs in reference namelist READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) IF(lwm) WRITE ( numond, namobs ) lk_diaobs = .FALSE. #if defined key_diaobs IF ( ln_diaobs ) lk_diaobs = .TRUE. #endif IF ( .NOT. lk_diaobs ) THEN IF(lwp) WRITE(numout,cform_war) IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' RETURN ENDIF IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia WRITE(numout,*) ' Logical switch for surface non-diatom logchl obs ln_slchlnon = ', ln_slchlnon WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia WRITE(numout,*) ' Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup IF (ln_grid_search_lookup) & WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight ENDIF !----------------------------------------------------------------------- ! Set up list of observation types to be used ! and the files associated with each type !----------------------------------------------------------------------- nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & & ln_sfco2, ln_spco2 /) ) IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN IF(lwp) WRITE(numout,cform_war) IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & & ' are set to .FALSE. so turning off calls to dia_obs' nwarn = nwarn + 1 lk_diaobs = .FALSE. RETURN ENDIF IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes IF ( nproftypes > 0 ) THEN ALLOCATE( cobstypesprof(nproftypes) ) ALLOCATE( ifilesprof(nproftypes) ) ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) jtype = 0 IF (ln_t3d .OR. ln_s3d) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'prof' clproffiles(jtype,:) = cn_profbfiles ENDIF IF (ln_vel3d) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'vel' clproffiles(jtype,:) = cn_velfbfiles ENDIF IF (ln_plchltot) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'plchltot' clproffiles(jtype,:) = cn_plchltotfbfiles ENDIF IF (ln_pchltot) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'pchltot' clproffiles(jtype,:) = cn_pchltotfbfiles ENDIF IF (ln_pno3) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'pno3' clproffiles(jtype,:) = cn_pno3fbfiles ENDIF IF (ln_psi4) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'psi4' clproffiles(jtype,:) = cn_psi4fbfiles ENDIF IF (ln_ppo4) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'ppo4' clproffiles(jtype,:) = cn_ppo4fbfiles ENDIF IF (ln_pdic) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'pdic' clproffiles(jtype,:) = cn_pdicfbfiles ENDIF IF (ln_palk) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'palk' clproffiles(jtype,:) = cn_palkfbfiles ENDIF IF (ln_pph) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'pph' clproffiles(jtype,:) = cn_pphfbfiles ENDIF IF (ln_po2) THEN jtype = jtype + 1 cobstypesprof(jtype) = 'po2' clproffiles(jtype,:) = cn_po2fbfiles ENDIF CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) ENDIF IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes IF ( nsurftypes > 0 ) THEN ALLOCATE( cobstypessurf(nsurftypes) ) ALLOCATE( ifilessurf(nsurftypes) ) ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) ALLOCATE(n2dintsurf(nsurftypes)) ALLOCATE(ravglamscl(nsurftypes)) ALLOCATE(ravgphiscl(nsurftypes)) ALLOCATE(lfpindegs(nsurftypes)) ALLOCATE(llnightav(nsurftypes)) jtype = 0 IF (ln_sla) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sla' clsurffiles(jtype,:) = cn_slafbfiles ENDIF IF (ln_sst) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sst' clsurffiles(jtype,:) = cn_sstfbfiles ENDIF IF (ln_sic) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sic' clsurffiles(jtype,:) = cn_sicfbfiles ENDIF IF (ln_sss) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sss' clsurffiles(jtype,:) = cn_sssfbfiles ENDIF IF (ln_slchltot) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchltot' clsurffiles(jtype,:) = cn_slchltotfbfiles ENDIF IF (ln_slchldia) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchldia' clsurffiles(jtype,:) = cn_slchldiafbfiles ENDIF IF (ln_slchlnon) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchlnon' clsurffiles(jtype,:) = cn_slchlnonfbfiles ENDIF IF (ln_slchldin) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchldin' clsurffiles(jtype,:) = cn_slchldinfbfiles ENDIF IF (ln_slchlmic) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchlmic' clsurffiles(jtype,:) = cn_slchlmicfbfiles ENDIF IF (ln_slchlnan) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchlnan' clsurffiles(jtype,:) = cn_slchlnanfbfiles ENDIF IF (ln_slchlpic) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slchlpic' clsurffiles(jtype,:) = cn_slchlpicfbfiles ENDIF IF (ln_schltot) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'schltot' clsurffiles(jtype,:) = cn_schltotfbfiles ENDIF IF (ln_slphytot) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slphytot' clsurffiles(jtype,:) = cn_slphytotfbfiles ENDIF IF (ln_slphydia) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slphydia' clsurffiles(jtype,:) = cn_slphydiafbfiles ENDIF IF (ln_slphynon) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'slphynon' clsurffiles(jtype,:) = cn_slphynonfbfiles ENDIF IF (ln_sspm) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sspm' clsurffiles(jtype,:) = cn_sspmfbfiles ENDIF IF (ln_sfco2) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'sfco2' clsurffiles(jtype,:) = cn_sfco2fbfiles ENDIF IF (ln_spco2) THEN jtype = jtype + 1 cobstypessurf(jtype) = 'spco2' clsurffiles(jtype,:) = cn_spco2fbfiles ENDIF CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) DO jtype = 1, nsurftypes IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN IF ( nn_2dint_sla == -1 ) THEN n2dint_type = nn_2dint_default ELSE n2dint_type = nn_2dint_sla ENDIF ztype_avglamscl = rn_sla_avglamscl ztype_avgphiscl = rn_sla_avgphiscl ltype_fp_indegs = ln_sla_fp_indegs ltype_night = .FALSE. ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN IF ( nn_2dint_sst == -1 ) THEN n2dint_type = nn_2dint_default ELSE n2dint_type = nn_2dint_sst ENDIF ztype_avglamscl = rn_sst_avglamscl ztype_avgphiscl = rn_sst_avgphiscl ltype_fp_indegs = ln_sst_fp_indegs ltype_night = ln_sstnight ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN IF ( nn_2dint_sic == -1 ) THEN n2dint_type = nn_2dint_default ELSE n2dint_type = nn_2dint_sic ENDIF ztype_avglamscl = rn_sic_avglamscl ztype_avgphiscl = rn_sic_avgphiscl ltype_fp_indegs = ln_sic_fp_indegs ltype_night = .FALSE. ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN IF ( nn_2dint_sss == -1 ) THEN n2dint_type = nn_2dint_default ELSE n2dint_type = nn_2dint_sss ENDIF ztype_avglamscl = rn_sss_avglamscl ztype_avgphiscl = rn_sss_avgphiscl ltype_fp_indegs = ln_sss_fp_indegs ltype_night = .FALSE. ELSE n2dint_type = nn_2dint_default ztype_avglamscl = rn_default_avglamscl ztype_avgphiscl = rn_default_avgphiscl ltype_fp_indegs = ln_default_fp_indegs ltype_night = .FALSE. ENDIF CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & & nn_2dint_default, n2dint_type, & & ztype_avglamscl, ztype_avgphiscl, & & ltype_fp_indegs, ltype_night, & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) END DO ENDIF IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' !----------------------------------------------------------------------- ! Obs operator parameter checking and initialisations !----------------------------------------------------------------------- IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) RETURN ENDIF IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN CALL ctl_stop(' Choice of vertical (1D) interpolation method', & & ' is not available') ENDIF IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & & ' is not available') ENDIF CALL obs_typ_init CALL mppmap_init CALL obs_grid_setup( ) !----------------------------------------------------------------------- ! Depending on switches read the various observation types !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN ALLOCATE(profdata(nproftypes)) ALLOCATE(profdataqc(nproftypes)) ALLOCATE(nvarsprof(nproftypes)) ALLOCATE(nextrprof(nproftypes)) DO jtype = 1, nproftypes IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN nvarsprof(jtype) = 2 nextrprof(jtype) = 1 ALLOCATE(llvar(nvarsprof(jtype))) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) llvar(1) = ln_t3d llvar(2) = ln_s3d zglam(:,:,1) = glamt(:,:) zglam(:,:,2) = glamt(:,:) zgphi(:,:,1) = gphit(:,:) zgphi(:,:,2) = gphit(:,:) zmask(:,:,:,1) = tmask(:,:,:) zmask(:,:,:,2) = tmask(:,:,:) ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN nvarsprof(jtype) = 2 nextrprof(jtype) = 2 ALLOCATE(llvar(nvarsprof(jtype))) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) llvar(1) = ln_vel3d llvar(2) = ln_vel3d zglam(:,:,1) = glamu(:,:) zglam(:,:,2) = glamv(:,:) zgphi(:,:,1) = gphiu(:,:) zgphi(:,:,2) = gphiv(:,:) zmask(:,:,:,1) = umask(:,:,:) zmask(:,:,:,2) = vmask(:,:,:) ELSE nvarsprof(jtype) = 1 nextrprof(jtype) = 0 ALLOCATE(llvar(nvarsprof(jtype))) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) llvar(1) = .TRUE. zglam(:,:,1) = glamt(:,:) zgphi(:,:,1) = gphit(:,:) zmask(:,:,:,1) = tmask(:,:,:) ENDIF !Read in profile or profile obs types CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & & clproffiles(jtype,1:ifilesprof(jtype)), & & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & & rn_dobsini, rn_dobsend, llvar, & & ln_ignmis, ln_s_at_t, .FALSE., & & kdailyavtypes = nn_profdavtypes ) DO jvar = 1, nvarsprof(jtype) CALL obs_prof_staend( profdata(jtype), jvar ) END DO CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & & llvar, & & jpi, jpj, jpk, & & zmask, zglam, zgphi, & & ln_nea, ln_bound_reject, & & kdailyavtypes = nn_profdavtypes ) DEALLOCATE( llvar ) CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zglam ) CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zgphi ) CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) END DO DEALLOCATE( ifilesprof, clproffiles ) ENDIF IF ( nsurftypes > 0 ) THEN ALLOCATE(surfdata(nsurftypes)) ALLOCATE(surfdataqc(nsurftypes)) ALLOCATE(nvarssurf(nsurftypes)) ALLOCATE(nextrsurf(nsurftypes)) DO jtype = 1, nsurftypes nvarssurf(jtype) = 1 nextrsurf(jtype) = 0 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 !Read in surface obs types CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & & clsurffiles(jtype,1:ifilessurf(jtype)), & & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) IF ( ln_altbias ) & & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) ENDIF IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN jnumsstbias = 0 DO jfile = 1, jpmaxnfiles IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & & jnumsstbias = jnumsstbias + 1 END DO IF ( jnumsstbias == 0 ) THEN CALL ctl_stop("ln_sstbias set but no bias files to read in") ENDIF CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) ENDIF END DO DEALLOCATE( ifilessurf, clsurffiles ) ENDIF END SUBROUTINE dia_obs_init SUBROUTINE dia_obs( kstp ) !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs *** !! !! ** Purpose : Call the observation operators on each time step !! !! ** Method : Call the observation operators on each time step to !! compute the model equivalent of the following data: !! - Profile data, currently T/S or U/V !! - Surface data, currently SST, SLA or sea-ice concentration. !! !! ** Action : !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 07-04 (G. Smith) Generalized surface operators !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles !! ! 15-08 (M. Martin) Combined surface/profile routines. !!---------------------------------------------------------------------- !! * Modules used USE phycst, ONLY : & ! Physical constants & rday USE oce, ONLY : & ! Ocean dynamics and tracers variables & tsn, & & un, & & vn, & & sshn #if defined key_lim3 USE ice, ONLY : & ! LIM3 Ice model variables & frld #endif #if defined key_lim2 USE ice_2, ONLY : & ! LIM2 Ice model variables & frld #endif #if defined key_cice USE sbc_oce, ONLY : fr_i ! ice fraction #endif #if defined key_hadocc USE trc, ONLY : & ! HadOCC variables & trn, & & HADOCC_CHL, & & HADOCC_FCO2, & & HADOCC_PCO2, & & HADOCC_FILL_FLT USE par_hadocc USE had_bgc_const, ONLY: c2n_p #elif defined key_medusa USE trc, ONLY : & ! MEDUSA variables & trn USE par_medusa USE sms_medusa, ONLY: & & xthetapn, & & xthetapd #if defined key_roam USE sms_medusa, ONLY: & & f2_pco2w, & & f2_fco2w, & & f3_pH #endif #elif defined key_fabm USE fabm USE par_fabm #endif #if defined key_spm USE par_spm, ONLY: & ! ERSEM/SPM sediments & jp_spm USE trc, ONLY : & & trn #endif IMPLICIT NONE !! * Arguments INTEGER, INTENT(IN) :: kstp ! Current timestep !! * Local declarations INTEGER :: idaystp ! Number of timesteps per day INTEGER :: jtype ! Data loop variable INTEGER :: jvar ! Variable number INTEGER :: ji, jj, jk ! Loop counters REAL(wp) :: tiny ! small number REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & & zprofvar ! Model values for variables in a prof ob REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & & zprofmask ! Mask associated with zprofvar REAL(wp), POINTER, DIMENSION(:,:) :: & & zsurfvar, & ! Model values equivalent to surface ob. & zsurfmask ! Mask associated with surface variable REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zglam, & ! Model longitudes for prof variables & zgphi ! Model latitudes for prof variables LOGICAL :: llog10 ! Perform log10 transform of variable IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'dia_obs : Call the observation operators', kstp WRITE(numout,*) '~~~~~~~' CALL FLUSH(numout) ENDIF idaystp = NINT( rday / rdt ) !----------------------------------------------------------------------- ! Call the profile and surface observation operators !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN DO jtype = 1, nproftypes ! Allocate local work arrays CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) ! Defaults which might change DO jvar = 1, profdataqc(jtype)%nvar zprofmask(:,:,:,jvar) = tmask(:,:,:) zglam(:,:,jvar) = glamt(:,:) zgphi(:,:,jvar) = gphit(:,:) END DO SELECT CASE ( TRIM(cobstypesprof(jtype)) ) CASE('prof') zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) CASE('vel') zprofvar(:,:,:,1) = un(:,:,:) zprofvar(:,:,:,2) = vn(:,:,:) zprofmask(:,:,:,1) = umask(:,:,:) zprofmask(:,:,:,2) = vmask(:,:,:) zglam(:,:,1) = glamu(:,:) zglam(:,:,2) = glamv(:,:) zgphi(:,:,1) = gphiu(:,:) zgphi(:,:,2) = gphiv(:,:) CASE('plchltot') #if defined key_hadocc ! Chlorophyll from HadOCC zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) #elif defined key_medusa ! Add non-diatom and diatom chlorophyll from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) #elif defined key_fabm ! Add all chlorophyll groups from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run plchltot observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif ! Take the log10 where we can, otherwise exclude tiny = 1.0e-20 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) ELSEWHERE zprofvar(:,:,:,:) = obfillflt zprofmask(:,:,:,:) = 0 END WHERE ! Mask out model below any excluded values, ! to avoid interpolation issues DO jvar = 1, profdataqc(jtype)%nvar DO jj = 1, jpj DO ji = 1, jpi depth_loop: DO jk = 1, jpk IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN zprofmask(ji,jj,jk:jpk,jvar) = 0 EXIT depth_loop ENDIF END DO depth_loop END DO END DO END DO CASE('pchltot') #if defined key_hadocc ! Chlorophyll from HadOCC zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) #elif defined key_medusa ! Add non-diatom and diatom chlorophyll from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) #elif defined key_fabm ! Add all chlorophyll groups from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & & trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run pchltot observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('pno3') #if defined key_hadocc ! Dissolved inorganic nitrogen from HadOCC zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) #elif defined key_medusa ! Dissolved inorganic nitrogen from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpdin) #elif defined key_fabm ! Nitrate from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n3n) #else CALL ctl_stop( ' Trying to run pno3 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('psi4') #if defined key_hadocc CALL ctl_stop( ' Trying to run psi4 observation operator', & & ' but HadOCC does not simulate silicate' ) #elif defined key_medusa ! Silicate from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpsil) #elif defined key_fabm ! Silicate from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n5s) #else CALL ctl_stop( ' Trying to run psi4 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('ppo4') #if defined key_hadocc CALL ctl_stop( ' Trying to run ppo4 observation operator', & & ' but HadOCC does not simulate phosphate' ) #elif defined key_medusa CALL ctl_stop( ' Trying to run ppo4 observation operator', & & ' but MEDUSA does not simulate phosphate' ) #elif defined key_fabm ! Phosphate from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n1p) #else CALL ctl_stop( ' Trying to run ppo4 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('pdic') #if defined key_hadocc ! Dissolved inorganic carbon from HadOCC zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) #elif defined key_medusa ! Dissolved inorganic carbon from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpdic) #elif defined key_fabm ! Dissolved inorganic carbon from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3c) #else CALL ctl_stop( ' Trying to run pdic observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('palk') #if defined key_hadocc ! Alkalinity from HadOCC zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) #elif defined key_medusa ! Alkalinity from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpalk) #elif defined key_fabm ! Alkalinity from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3a) #else CALL ctl_stop( ' Trying to run palk observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('pph') #if defined key_hadocc CALL ctl_stop( ' Trying to run pph observation operator', & & ' but HadOCC has no pH diagnostic defined' ) #elif defined key_medusa && defined key_roam ! pH from MEDUSA zprofvar(:,:,:,1) = f3_pH(:,:,:) #elif defined key_fabm ! pH from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3ph) #else CALL ctl_stop( ' Trying to run pph observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('po2') #if defined key_hadocc CALL ctl_stop( ' Trying to run po2 observation operator', & & ' but HadOCC does not simulate oxygen' ) #elif defined key_medusa ! Oxygen from MEDUSA zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) #elif defined key_fabm ! Oxygen from ERSEM zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o2o) #else CALL ctl_stop( ' Trying to run po2 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE DEFAULT CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) END SELECT DO jvar = 1, profdataqc(jtype)%nvar CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & & nit000, idaystp, jvar, & & zprofvar(:,:,:,jvar), & & fsdept(:,:,:), fsdepw(:,:,:), & & zprofmask(:,:,:,jvar), & & zglam(:,:,jvar), zgphi(:,:,jvar), & & nn_1dint, nn_2dint_default, & & kdailyavtypes = nn_profdavtypes ) END DO CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) END DO ENDIF IF ( nsurftypes > 0 ) THEN !Allocate local work arrays CALL wrk_alloc( jpi, jpj, zsurfvar ) CALL wrk_alloc( jpi, jpj, zsurfmask ) DO jtype = 1, nsurftypes !Defaults which might be changed zsurfmask(:,:) = tmask(:,:,1) llog10 = .FALSE. SELECT CASE ( TRIM(cobstypessurf(jtype)) ) CASE('sst') zsurfvar(:,:) = tsn(:,:,1,jp_tem) CASE('sla') zsurfvar(:,:) = sshn(:,:) CASE('sss') zsurfvar(:,:) = tsn(:,:,1,jp_sal) CASE('sic') IF ( kstp == 0 ) THEN IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & & 'time-step but some obs are valid then.' ) WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & & ' sea-ice obs will be missed' ENDIF surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & & surfdataqc(jtype)%nsstp(1) CYCLE ELSE #if defined key_cice zsurfvar(:,:) = fr_i(:,:) #elif defined key_lim2 || defined key_lim3 zsurfvar(:,:) = 1._wp - frld(:,:) #else CALL ctl_stop( ' Trying to run sea-ice observation operator', & & ' but no sea-ice model appears to have been defined' ) #endif ENDIF CASE('slchltot') #if defined key_hadocc ! Surface chlorophyll from HadOCC zsurfvar(:,:) = HADOCC_CHL(:,:,1) #elif defined key_medusa ! Add non-diatom and diatom surface chlorophyll from MEDUSA zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) #elif defined key_fabm ! Add all surface chlorophyll groups from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl2) + & & trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run slchltot observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchldia') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchldia observation operator', & & ' but HadOCC does not explicitly simulate diatoms' ) #elif defined key_medusa ! Diatom surface chlorophyll from MEDUSA zsurfvar(:,:) = trn(:,:,1,jpchd) #elif defined key_fabm ! Diatom surface chlorophyll from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) #else CALL ctl_stop( ' Trying to run slchldia observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchlnon') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchlnon observation operator', & & ' but HadOCC does not explicitly simulate non-diatoms' ) #elif defined key_medusa ! Non-diatom surface chlorophyll from MEDUSA zsurfvar(:,:) = trn(:,:,1,jpchn) #elif defined key_fabm ! Add all non-diatom surface chlorophyll groups from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl2) + & & trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run slchlnon observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchldin') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchldin observation operator', & & ' but HadOCC does not explicitly simulate dinoflagellates' ) #elif defined key_medusa CALL ctl_stop( ' Trying to run slchldin observation operator', & & ' but MEDUSA does not explicitly simulate dinoflagellates' ) #elif defined key_fabm ! Dinoflagellate surface chlorophyll from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run slchldin observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchlmic') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchlmic observation operator', & & ' but HadOCC does not explicitly simulate microphytoplankton' ) #elif defined key_medusa CALL ctl_stop( ' Trying to run slchlmic observation operator', & & ' but MEDUSA does not explicitly simulate microphytoplankton' ) #elif defined key_fabm ! Add diatom and dinoflagellate surface chlorophyll from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run slchlmic observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchlnan') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchlnan observation operator', & & ' but HadOCC does not explicitly simulate nanophytoplankton' ) #elif defined key_medusa CALL ctl_stop( ' Trying to run slchlnan observation operator', & & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) #elif defined key_fabm ! Nanophytoplankton surface chlorophyll from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl2) #else CALL ctl_stop( ' Trying to run slchlnan observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slchlpic') #if defined key_hadocc CALL ctl_stop( ' Trying to run slchlpic observation operator', & & ' but HadOCC does not explicitly simulate picophytoplankton' ) #elif defined key_medusa CALL ctl_stop( ' Trying to run slchlpic observation operator', & & ' but MEDUSA does not explicitly simulate picophytoplankton' ) #elif defined key_fabm ! Picophytoplankton surface chlorophyll from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl3) #else CALL ctl_stop( ' Trying to run slchlpic observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('schltot') #if defined key_hadocc ! Surface chlorophyll from HadOCC zsurfvar(:,:) = HADOCC_CHL(:,:,1) #elif defined key_medusa ! Add non-diatom and diatom surface chlorophyll from MEDUSA zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) #elif defined key_fabm ! Add all surface chlorophyll groups from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl2) + & & trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4) #else CALL ctl_stop( ' Trying to run schltot observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('slphytot') #if defined key_hadocc ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p #elif defined key_medusa ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA ! multiplied by C:N ratio for each zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) #elif defined key_fabm ! Add all surface phytoplankton carbon groups from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_p1c) + trn(:,:,1,jp_fabm_p2c) + & & trn(:,:,1,jp_fabm_p3c) + trn(:,:,1,jp_fabm_p4c) #else CALL ctl_stop( ' Trying to run slphytot observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slphydia') #if defined key_hadocc CALL ctl_stop( ' Trying to run slphydia observation operator', & & ' but HadOCC does not explicitly simulate diatoms' ) #elif defined key_medusa ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd #elif defined key_fabm ! Diatom surface phytoplankton carbon from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_p1c) #else CALL ctl_stop( ' Trying to run slphydia observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('slphynon') #if defined key_hadocc CALL ctl_stop( ' Trying to run slphynon observation operator', & & ' but HadOCC does not explicitly simulate non-diatoms' ) #elif defined key_medusa ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn #elif defined key_fabm ! Add all non-diatom surface phytoplankton carbon groups from ERSEM zsurfvar(:,:) = trn(:,:,1,jp_fabm_p2c) + & & trn(:,:,1,jp_fabm_p3c) + trn(:,:,1,jp_fabm_p4c) #else CALL ctl_stop( ' Trying to run slphynon observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif llog10 = .TRUE. CASE('sspm') #if defined key_spm zsurfvar(:,:) = 0.0 DO jn = 1, jp_spm zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes END DO #else CALL ctl_stop( ' Trying to run sspm observation operator', & & ' but no spm model appears to have been defined' ) #endif CASE('sfco2') #if defined key_hadocc zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) ENDIF #elif defined key_medusa && defined key_roam zsurfvar(:,:) = f2_fco2w(:,:) #elif defined key_fabm ! First, get pCO2 from FABM pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) zsurfvar(:,:) = pco2_3d(:,:,1) ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems ! and data reduction routines, Deep-Sea Research II, 56: 512-522. ! and ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, ! Marine Chemistry, 2: 203-215. ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so ! not explicitly included - atmospheric pressure is not necessarily available so this is ! the best assumption. ! Further, the (1-xCO2)^2 term has been neglected. This is common practice ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) #else CALL ctl_stop( ' Trying to run sfco2 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('spco2') #if defined key_hadocc zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) ENDIF #elif defined key_medusa && defined key_roam zsurfvar(:,:) = f2_pco2w(:,:) #elif defined key_fabm pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) zsurfvar(:,:) = pco2_3d(:,:,1) #else CALL ctl_stop( ' Trying to run spco2 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE DEFAULT CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) END SELECT IF ( llog10 ) THEN ! Take the log10 where we can, otherwise exclude tiny = 1.0e-20 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) zsurfvar(:,:) = LOG10(zsurfvar(:,:)) ELSEWHERE zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 END WHERE ENDIF CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & & nit000, idaystp, zsurfvar, zsurfmask, & & n2dintsurf(jtype), llnightav(jtype), & & ravglamscl(jtype), ravgphiscl(jtype), & & lfpindegs(jtype) ) END DO CALL wrk_dealloc( jpi, jpj, zsurfvar ) CALL wrk_dealloc( jpi, jpj, zsurfmask ) ENDIF END SUBROUTINE dia_obs SUBROUTINE dia_obs_wri !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_wri *** !! !! ** Purpose : Call observation diagnostic output routines !! !! ** Method : Call observation diagnostic output routines !! !! ** Action : !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles !! ! 15-08 (M. Martin) Combined writing for prof and surf types !!---------------------------------------------------------------------- !! * Modules used USE obs_rot_vel ! Rotation of velocities IMPLICIT NONE !! * Local declarations INTEGER :: jtype ! Data set loop variable INTEGER :: jo, jvar, jk REAL(wp), DIMENSION(:), ALLOCATABLE :: & & zu, & & zv !----------------------------------------------------------------------- ! Depending on switches call various observation output routines !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN DO jtype = 1, nproftypes IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN ! For velocity data, rotate the model velocities to N/S, E/W ! using the compressed data structure. ALLOCATE( & & zu(profdataqc(jtype)%nvprot(1)), & & zv(profdataqc(jtype)%nvprot(2)) & & ) CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) DO jo = 1, profdataqc(jtype)%nprof DO jvar = 1, 2 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) IF ( jvar == 1 ) THEN profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) ELSE profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) ENDIF END DO END DO END DO DEALLOCATE( zu ) DEALLOCATE( zv ) END IF CALL obs_prof_decompress( profdataqc(jtype), & & profdata(jtype), .TRUE., numout ) CALL obs_wri_prof( profdata(jtype) ) END DO ENDIF IF ( nsurftypes > 0 ) THEN DO jtype = 1, nsurftypes CALL obs_surf_decompress( surfdataqc(jtype), & & surfdata(jtype), .TRUE., numout ) CALL obs_wri_surf( surfdata(jtype) ) END DO ENDIF END SUBROUTINE dia_obs_wri SUBROUTINE dia_obs_dealloc IMPLICIT NONE !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_dealloc *** !! !! ** Purpose : To deallocate data to enable the obs_oper online loop. !! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri !! !! ** Method : Clean up various arrays left behind by the obs_oper. !! !! ** Action : !! !!---------------------------------------------------------------------- ! obs_grid deallocation CALL obs_grid_deallocate ! diaobs deallocation IF ( nproftypes > 0 ) & & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) IF ( nsurftypes > 0 ) & & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) END SUBROUTINE dia_obs_dealloc SUBROUTINE ini_date( ddobsini ) !!---------------------------------------------------------------------- !! *** ROUTINE ini_date *** !! !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 06-10 (G. Smith) Calculates initial date the same as method for final date !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 !!---------------------------------------------------------------------- USE phycst, ONLY : & ! Physical constants & rday USE dom_oce, ONLY : & ! Ocean space and time domain variables & rdt IMPLICIT NONE !! * Arguments REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS !! * Local declarations INTEGER :: iyea ! date - (year, month, day, hour, minute) INTEGER :: imon INTEGER :: iday INTEGER :: ihou INTEGER :: imin INTEGER :: imday ! Number of days in month. INTEGER, DIMENSION(12) :: & & imonth_len ! Length in days of the months of the current year REAL(wp) :: zdayfrc ! Fraction of day !---------------------------------------------------------------------- ! Initial date initialization (year, month, day, hour, minute) ! (This assumes that the initial date is for 00z)) !---------------------------------------------------------------------- iyea = ndate0 / 10000 imon = ( ndate0 - iyea * 10000 ) / 100 iday = ndate0 - iyea * 10000 - imon * 100 ihou = 0 imin = 0 !---------------------------------------------------------------------- ! Compute number of days + number of hours + min since initial time !---------------------------------------------------------------------- iday = iday + ( nit000 -1 ) * rdt / rday zdayfrc = ( nit000 -1 ) * rdt / rday zdayfrc = zdayfrc - aint(zdayfrc) ihou = int( zdayfrc * 24 ) imin = int( (zdayfrc * 24 - ihou) * 60 ) !----------------------------------------------------------------------- ! Convert number of days (iday) into a real date !---------------------------------------------------------------------- CALL calc_month_len( iyea, imonth_len ) DO WHILE ( iday > imonth_len(imon) ) iday = iday - imonth_len(imon) imon = imon + 1 IF ( imon > 12 ) THEN imon = 1 iyea = iyea + 1 CALL calc_month_len( iyea, imonth_len ) ! update month lengths ENDIF END DO !---------------------------------------------------------------------- ! Convert it into YYYYMMDD.HHMMSS format. !---------------------------------------------------------------------- ddobsini = iyea * 10000_dp + imon * 100_dp + & & iday + ihou * 0.01_dp + imin * 0.0001_dp END SUBROUTINE ini_date SUBROUTINE fin_date( ddobsfin ) !!---------------------------------------------------------------------- !! *** ROUTINE fin_date *** !! !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format !! !! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format !! !! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 !!---------------------------------------------------------------------- USE phycst, ONLY : & ! Physical constants & rday USE dom_oce, ONLY : & ! Ocean space and time domain variables & rdt IMPLICIT NONE !! * Arguments REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS !! * Local declarations INTEGER :: iyea ! date - (year, month, day, hour, minute) INTEGER :: imon INTEGER :: iday INTEGER :: ihou INTEGER :: imin INTEGER :: imday ! Number of days in month. INTEGER, DIMENSION(12) :: & & imonth_len ! Length in days of the months of the current year REAL(wp) :: zdayfrc ! Fraction of day !----------------------------------------------------------------------- ! Initial date initialization (year, month, day, hour, minute) ! (This assumes that the initial date is for 00z) !----------------------------------------------------------------------- iyea = ndate0 / 10000 imon = ( ndate0 - iyea * 10000 ) / 100 iday = ndate0 - iyea * 10000 - imon * 100 ihou = 0 imin = 0 !----------------------------------------------------------------------- ! Compute number of days + number of hours + min since initial time !----------------------------------------------------------------------- iday = iday + nitend * rdt / rday zdayfrc = nitend * rdt / rday zdayfrc = zdayfrc - AINT( zdayfrc ) ihou = INT( zdayfrc * 24 ) imin = INT( ( zdayfrc * 24 - ihou ) * 60 ) !----------------------------------------------------------------------- ! Convert number of days (iday) into a real date !---------------------------------------------------------------------- CALL calc_month_len( iyea, imonth_len ) DO WHILE ( iday > imonth_len(imon) ) iday = iday - imonth_len(imon) imon = imon + 1 IF ( imon > 12 ) THEN imon = 1 iyea = iyea + 1 CALL calc_month_len( iyea, imonth_len ) ! update month lengths ENDIF END DO !----------------------------------------------------------------------- ! Convert it into YYYYMMDD.HHMMSS format !----------------------------------------------------------------------- ddobsfin = iyea * 10000_dp + imon * 100_dp + iday & & + ihou * 0.01_dp + imin * 0.0001_dp END SUBROUTINE fin_date SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) INTEGER, INTENT(IN) :: ntypes ! Total number of obs types INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & & ifiles ! Out number of files for each type CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & & cobstypes ! List of obs types CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & & cfiles ! List of files for all types !Local variables INTEGER :: jfile INTEGER :: jtype DO jtype = 1, ntypes ifiles(jtype) = 0 DO jfile = 1, jpmaxnfiles IF ( trim(cfiles(jtype,jfile)) /= '' ) & ifiles(jtype) = ifiles(jtype) + 1 END DO IF ( ifiles(jtype) == 0 ) THEN CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// & & ' set to true but no files available to read' ) ENDIF IF(lwp) THEN WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' DO jfile = 1, ifiles(jtype) WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) END DO ENDIF END DO END SUBROUTINE obs_settypefiles SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & & n2dint_default, n2dint_type, & & ravglamscl_type, ravgphiscl_type, & & lfp_indegs_type, lavnight_type, & & n2dint, ravglamscl, ravgphiscl, & & lfpindegs, lavnight ) INTEGER, INTENT(IN) :: ntypes ! Total number of obs types INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type REAL(wp), INTENT(IN) :: & & ravglamscl_type, & !E/W diameter of obs footprint for this type & ravgphiscl_type !N/S diameter of obs footprint for this type LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average CHARACTER(len=8), INTENT(IN) :: ctypein INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & & n2dint REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & & ravglamscl, ravgphiscl LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & & lfpindegs, lavnight lavnight(jtype) = lavnight_type IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN n2dint(jtype) = n2dint_type ELSE IF ( n2dint_type == -1 ) THEN n2dint(jtype) = n2dint_default ELSE CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & & ' is not available') ENDIF ! For averaging observation footprints set options for size of footprint IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN IF ( ravglamscl_type > 0._wp ) THEN ravglamscl(jtype) = ravglamscl_type ELSE CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) ENDIF IF ( ravgphiscl_type > 0._wp ) THEN ravgphiscl(jtype) = ravgphiscl_type ELSE CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) ENDIF lfpindegs(jtype) = lfp_indegs_type ENDIF ! Write out info IF(lwp) THEN IF ( n2dint(jtype) <= 4 ) THEN WRITE(numout,*) ' '//TRIM(ctypein)// & & ' model counterparts will be interpolated horizontally' ELSE IF ( n2dint(jtype) <= 6 ) THEN WRITE(numout,*) ' '//TRIM(ctypein)// & & ' model counterparts will be averaged horizontally' WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) IF ( lfpindegs(jtype) ) THEN WRITE(numout,*) ' '//' (in degrees)' ELSE WRITE(numout,*) ' '//' (in metres)' ENDIF ENDIF ENDIF END SUBROUTINE obs_setinterpopts END MODULE diaobs