Changeset 15670
- Timestamp:
- 2022-01-25T15:20:24+01:00 (3 years ago)
- Location:
- branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 9 added
- 20 deleted
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r10728 r15670 76 76 LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments 77 77 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 78 LOGICAL, PUBLIC :: ln_ssh_hs_cons = .FALSE. !: Conserve heat and salt when adding SSH increment 78 79 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 79 80 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment … … 88 89 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 89 90 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 90 #if defined key_asminc91 91 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment 92 #endif93 92 ! !!! time steps relative to the cycle interval [0,nitend-nit000-1] 94 93 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term … … 173 172 & ln_pno3inc, ln_psi4inc, ln_pdicinc, ln_palkinc, & 174 173 & ln_pphinc, ln_po2inc, ln_ppo4inc, & 175 & ln_asmdin, ln_asmiau, 174 & ln_asmdin, ln_asmiau, ln_ssh_hs_cons, & 176 175 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 177 176 & ln_salfix, salfixmin, nn_divdmp, nitavgbkg, & … … 193 192 ln_asmiau = .TRUE. 194 193 ln_salfix = .FALSE. 194 ln_ssh_hs_cons = .FALSE. 195 195 ln_temnofreeze = .FALSE. 196 196 salfixmin = -9999 … … 222 222 WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc 223 223 WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc 224 WRITE(numout,*) ' Logical switch for conserving heat/salt when applying SSH increments ln_ssh_hs_cons = ', ln_ssh_hs_cons 224 225 WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc 225 226 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90
r10390 r15670 15 15 USE diurnal_bulk 16 16 USE cool_skin 17 #if defined key_fabm 18 USE par_fabm 19 USE fabm, ONLY: fabm_get_bulk_diagnostic_data 20 #endif 17 21 18 22 … … 109 113 CALL iom_put( "voce_op" , vn ) ! j-current 110 114 !CALL iom_put( "woce_op" , wn ) ! k-current 115 #if defined key_spm 116 cltra = TRIM(ctrc3d(5))//"_op" 117 zw3d(:,:,:) = trc3d(:,:,:,5)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) ! Visibility 118 CALL iom_put( cltra, zw3d ) 119 #endif 120 #if defined key_fabm 121 zw3d(:,:,:) = (1.7/fabm_get_bulk_diagnostic_data(model, jp_fabm_xeps))*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) ! hourly visibility 122 CALL iom_put( "Visib_op" , zw3d(:,:,:) ) ! hourly visibility 123 #endif 111 124 CALL calc_max_cur(zwu,zwv,zwz,zmdi) 112 125 CALL iom_put( "maxu" , zwu ) ! max u current -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r8058 r15670 6 6 !!====================================================================== 7 7 8 !!----------------------------------------------------------------------9 !! 'key_diaobs' : Switch on the observation diagnostic computation10 8 !!---------------------------------------------------------------------- 11 9 !! dia_obs_init : Reading and prepare observations … … 15 13 !! fin_date : Compute the final date YYYYMMDD.HHMMSS 16 14 !!---------------------------------------------------------------------- 17 !! * Modules used 15 !! * Modules used 18 16 USE wrk_nemo ! Memory Allocation 19 17 USE par_kind ! Precision variables … … 21 19 USE par_oce 22 20 USE dom_oce ! Ocean space and time domain variables 23 USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch 24 USE obs_read_prof ! Reading and allocation of observations (Coriolis) 25 USE obs_read_sla ! Reading and allocation of SLA observations 26 USE obs_read_sst ! Reading and allocation of SST observations 21 USE obs_read_prof ! Reading and allocation of profile obs 22 USE obs_read_surf ! Reading and allocation of surface obs 27 23 USE obs_readmdt ! Reading and allocation of MDT for SLA. 28 USE obs_read_seaice ! Reading and allocation of Sea Ice observations29 USE obs_read_vel ! Reading and allocation of velocity component observations30 24 USE obs_prep ! Preparation of obs. (grid search etc). 31 25 USE obs_oper ! Observation operators … … 33 27 USE obs_grid ! Grid searching 34 28 USE obs_read_altbias ! Bias treatment for altimeter 29 USE obs_sstbias ! Bias correction routine for SST 35 30 USE obs_profiles_def ! Profile data definitions 36 USE obs_profiles ! Profile data storage37 31 USE obs_surf_def ! Surface data definitions 38 USE obs_sla ! SLA data storage39 USE obs_sst ! SST data storage40 USE obs_seaice ! Sea Ice data storage41 32 USE obs_types ! Definitions for observation types 42 33 USE mpp_map ! MPP mapping 43 34 USE lib_mpp ! For ctl_warn/stop 35 USE tradmp ! For climatological temperature & salinity 44 36 45 37 IMPLICIT NONE … … 52 44 & dia_obs_dealloc ! Deallocate dia_obs data 53 45 54 !! * Shared Module variables55 LOGICAL, PUBLIC, PARAMETER :: &56 #if defined key_diaobs57 & lk_diaobs = .TRUE. !: Logical switch for observation diangostics58 #else59 & lk_diaobs = .FALSE. !: Logical switch for observation diangostics60 #endif61 62 46 !! * Module variables 63 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles64 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles65 LOGICAL , PUBLIC :: ln_ena !: Logical switch for the ENACT data set66 LOGICAL , PUBLIC :: ln_cor !: Logical switch for the Coriolis data set67 LOGICAL , PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles68 LOGICAL , PUBLIC :: ln_sla !: Logical switch for sea level anomalies69 LOGICAL , PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files70 LOGICAL , PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files71 LOGICAL , PUBLIC :: ln_sst !: Logical switch for sea surface temperature72 LOGICAL , PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature73 LOGICAL , PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data74 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files 75 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration76 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations77 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data78 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data79 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data80 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data81 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files82 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height83 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity84 LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations85 LOGICAL, PUBLIC :: ln_nea !: Remove observations near land86 LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias87 LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files 88 LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations 89 90 REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS91 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS92 93 INTEGER , PUBLIC :: n1dint !: Vertical interpolation method94 INTEGER , PUBLIC :: n2dint !: Horizontal interpolation method95 47 LOGICAL, PUBLIC :: & 48 & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. 49 LOGICAL :: ln_diaobs !: Logical switch for the obs operator 50 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 51 LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 52 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 54 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 55 LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 56 LOGICAL :: ln_output_clim !: Logical switch for interpolating and writing T/S climatology 57 LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 58 59 REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 60 REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 61 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint 62 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint 63 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint 64 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint 65 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint 66 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint 67 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint 68 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint 69 REAL(wp), PUBLIC :: & 70 & MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 71 72 73 INTEGER :: nn_1dint !: Vertical interpolation method 74 INTEGER :: nn_2dint_default !: Default horizontal interpolation method 75 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (-1 = default) 76 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (-1 = default) 77 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (-1 = default) 78 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method (-1 = default) 79 96 80 INTEGER, DIMENSION(imaxavtypes) :: & 97 & endailyavtypes !: ENACT data types which are daily average 98 99 INTEGER, PARAMETER :: MaxNumFiles = 1000 100 LOGICAL, DIMENSION(MaxNumFiles) :: & 101 & ln_profb_ena, & !: Is the feedback files from ENACT data ? 102 ! !: If so use endailyavtypes 103 & ln_profb_enatim !: Change tim for 820 enact data set. 104 105 LOGICAL, DIMENSION(MaxNumFiles) :: & 106 & ln_velfb_av !: Is the velocity feedback files daily average? 81 & nn_profdavtypes !: Profile data types representing a daily average 82 INTEGER :: nproftypes !: Number of profile obs types 83 INTEGER :: nsurftypes !: Number of surface obs types 84 INTEGER, DIMENSION(:), ALLOCATABLE :: & 85 & nvarsprof, & !: Number of profile variables 86 & nvarssurf !: Number of surface variables 87 INTEGER, DIMENSION(:), ALLOCATABLE :: & 88 & nextrprof, & !: Number of profile extra variables 89 & nextrsurf !: Number of surface extra variables 90 INTEGER, DIMENSION(:), ALLOCATABLE :: & 91 & n2dintsurf !: Interpolation option for surface variables 92 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 93 & ravglamscl, & !: E/W diameter of averaging footprint for surface variables 94 & ravgphiscl !: N/S diameter of averaging footprint for surface variables 107 95 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 108 & ld_enact !: Profile data is ENACT so use endailyavtypes 109 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 110 & ld_velav !: Velocity data is daily averaged 111 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 112 & ld_sstnight !: SST observation corresponds to night mean 96 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 97 & llnightav !: Logical for calculating night-time averages 98 99 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 100 & surfdata, & !: Initial surface data 101 & surfdataqc !: Surface data after quality control 102 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 103 & profdata, & !: Initial profile data 104 & profdataqc !: Profile data after quality control 105 106 CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 107 & cobstypesprof, & !: Profile obs types 108 & cobstypessurf !: Surface obs types 113 109 114 110 !!---------------------------------------------------------------------- … … 118 114 !!---------------------------------------------------------------------- 119 115 116 !! * Substitutions 117 # include "domzgr_substitute.h90" 120 118 CONTAINS 121 119 … … 135 133 !! ! 06-10 (A. Weaver) Cleaning and add controls 136 134 !! ! 07-03 (K. Mogensen) General handling of profiles 135 !! ! 14-08 (J.While) Incorporated SST bias correction 136 !! ! 15-02 (M. Martin) Simplification of namelist and code 137 137 !!---------------------------------------------------------------------- 138 138 … … 140 140 141 141 !! * Local declarations 142 CHARACTER(len=128) :: enactfiles(MaxNumFiles) 143 CHARACTER(len=128) :: coriofiles(MaxNumFiles) 144 CHARACTER(len=128) :: profbfiles(MaxNumFiles) 145 CHARACTER(len=128) :: sstfiles(MaxNumFiles) 146 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 147 CHARACTER(len=128) :: slafilesact(MaxNumFiles) 148 CHARACTER(len=128) :: slafilespas(MaxNumFiles) 149 CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 150 CHARACTER(len=128) :: seaicefiles(MaxNumFiles) 151 CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 152 CHARACTER(len=128) :: veladcpfiles(MaxNumFiles) 153 CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 154 CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 155 CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 156 CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 157 CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 158 CHARACTER(LEN=128) :: reysstname 159 CHARACTER(LEN=12) :: reysstfmt 160 CHARACTER(LEN=128) :: bias_file 161 CHARACTER(LEN=20) :: datestr=" ", timestr=" " 162 NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, & 163 & ln_sla, ln_sladt, ln_slafb, & 164 & ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, & 165 & enactfiles, coriofiles, profbfiles, & 166 & slafilesact, slafilespas, slafbfiles, & 167 & sstfiles, sstfbfiles, & 168 & ln_seaice, seaicefiles, & 169 & dobsini, dobsend, n1dint, n2dint, & 170 & nmsshc, mdtcorr, mdtcutoff, & 171 & ln_reysst, ln_ghrsst, reysstname, reysstfmt, & 172 & ln_sstnight, & 173 & ln_grid_search_lookup, & 174 & grid_search_file, grid_search_res, & 175 & ln_grid_global, bias_file, ln_altbias, & 176 & endailyavtypes, ln_s_at_t, ln_profb_ena, & 177 & ln_vel3d, ln_velavcur, velavcurfiles, & 178 & ln_velhrcur, velhrcurfiles, & 179 & ln_velavadcp, velavadcpfiles, & 180 & ln_velhradcp, velhradcpfiles, & 181 & ln_velfb, velfbfiles, ln_velfb_av, & 182 & ln_profb_enatim, ln_ignmis, ln_cl4 183 184 INTEGER :: jprofset 185 INTEGER :: jveloset 186 INTEGER :: jvar 187 INTEGER :: jnumenact 188 INTEGER :: jnumcorio 189 INTEGER :: jnumprofb 190 INTEGER :: jnumslaact 191 INTEGER :: jnumslapas 192 INTEGER :: jnumslafb 193 INTEGER :: jnumsst 194 INTEGER :: jnumsstfb 195 INTEGER :: jnumseaice 196 INTEGER :: jnumvelavcur 197 INTEGER :: jnumvelhrcur 198 INTEGER :: jnumvelavadcp 199 INTEGER :: jnumvelhradcp 200 INTEGER :: jnumvelfb 201 INTEGER :: ji 202 INTEGER :: jset 203 INTEGER :: ios ! Local integer output status for namelist read 204 LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 142 INTEGER, PARAMETER :: & 143 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 144 INTEGER, DIMENSION(:), ALLOCATABLE :: & 145 & ifilesprof, & ! Number of profile files 146 & ifilessurf ! Number of surface files 147 INTEGER :: ios ! Local integer output status for namelist read 148 INTEGER :: jtype ! Counter for obs types 149 INTEGER :: jvar ! Counter for variables 150 INTEGER :: jfile ! Counter for files 151 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 152 INTEGER :: n2dint_type ! Local version of nn_2dint* 153 154 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 155 & cn_profbfiles, & ! T/S profile input filenames 156 & cn_sstfbfiles, & ! Sea surface temperature input filenames 157 & cn_slafbfiles, & ! Sea level anomaly input filenames 158 & cn_sicfbfiles, & ! Seaice concentration input filenames 159 & cn_velfbfiles, & ! Velocity profile input filenames 160 & cn_sssfbfiles, & ! Sea surface salinity input filenames 161 & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames 162 & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames 163 & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames 164 & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames 165 & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 166 & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames 167 & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames 168 & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames 169 & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames 170 & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames 171 & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 172 & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames 173 & cn_skd490fbfiles, & ! Surface Kd490 input filenames 174 & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames 175 & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames 176 & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 177 & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames 178 & cn_pno3fbfiles, & ! Profile nitrate input filenames 179 & cn_psi4fbfiles, & ! Profile silicate input filenames 180 & cn_ppo4fbfiles, & ! Profile phosphate input filenames 181 & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames 182 & cn_palkfbfiles, & ! Profile alkalinity input filenames 183 & cn_pphfbfiles, & ! Profile pH input filenames 184 & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames 185 & cn_sstbiasfiles ! SST bias input filenames 186 187 CHARACTER(LEN=128) :: & 188 & cn_altbiasfile ! Altimeter bias input filename 189 190 191 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 192 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 193 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 194 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 195 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 196 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 197 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 198 LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs 199 LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs 200 LOGICAL :: ln_slchlnon ! Logical switch for surface non-diatom log10(chlorophyll) obs 201 LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs 202 LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 203 LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs 204 LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs 205 LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs 206 LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs 207 LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs 208 LOGICAL :: ln_slphynon ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 209 LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs 210 LOGICAL :: ln_skd490 ! Logical switch for surface Kd490 211 LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs 212 LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs 213 LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs 214 LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs 215 LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs 216 LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs 217 LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs 218 LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs 219 LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs 220 LOGICAL :: ln_pph ! Logical switch for profile pH obs 221 LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs 222 LOGICAL :: ln_nea ! Logical switch to remove obs near land 223 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 224 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 225 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 226 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 227 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 228 229 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 230 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 231 232 REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 233 REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 234 235 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 236 & clproffiles, & ! Profile filenames 237 & clsurffiles ! Surface filenames 238 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars ! Expected variable names 239 240 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read 241 LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 242 LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) 243 LOGICAL :: ltype_clim ! Local version of ln_output_clim 244 245 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 246 & zglam ! Model longitudes for profile variables 247 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 248 & zgphi ! Model latitudes for profile variables 249 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 250 & zmask ! Model land/sea mask associated with variables 251 252 253 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 254 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 255 & ln_slchltot, ln_slchldia, ln_slchlnon, & 256 & ln_slchldin, ln_slchlmic, ln_slchlnan, & 257 & ln_slchlpic, ln_schltot, & 258 & ln_slphytot, ln_slphydia, ln_slphynon, & 259 & ln_sspm, ln_sfco2, ln_spco2, & 260 & ln_skd490, & 261 & ln_plchltot, ln_pchltot, ln_pno3, & 262 & ln_psi4, ln_ppo4, ln_pdic, & 263 & ln_palk, ln_pph, ln_po2, & 264 & ln_altbias, ln_sstbias, ln_nea, & 265 & ln_grid_global, ln_grid_search_lookup, & 266 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 267 & ln_sstnight, ln_output_clim, & 268 & ln_time_mean_sla_bkg, ln_default_fp_indegs, & 269 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 270 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 271 & cn_profbfiles, cn_slafbfiles, & 272 & cn_sstfbfiles, cn_sicfbfiles, & 273 & cn_velfbfiles, cn_sssfbfiles, & 274 & cn_slchltotfbfiles, cn_slchldiafbfiles, & 275 & cn_slchlnonfbfiles, cn_slchldinfbfiles, & 276 & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & 277 & cn_slchlpicfbfiles, cn_schltotfbfiles, & 278 & cn_slphytotfbfiles, cn_slphydiafbfiles, & 279 & cn_slphynonfbfiles, cn_sspmfbfiles, & 280 & cn_skd490fbfiles, & 281 & cn_sfco2fbfiles, cn_spco2fbfiles, & 282 & cn_plchltotfbfiles, cn_pchltotfbfiles, & 283 & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 284 & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & 285 & cn_po2fbfiles, & 286 & cn_sstbiasfiles, cn_altbiasfile, & 287 & cn_gridsearchfile, rn_gridsearchres, & 288 & rn_dobsini, rn_dobsend, & 289 & rn_default_avglamscl, rn_default_avgphiscl, & 290 & rn_sla_avglamscl, rn_sla_avgphiscl, & 291 & rn_sst_avglamscl, rn_sst_avgphiscl, & 292 & rn_sss_avglamscl, rn_sss_avgphiscl, & 293 & rn_sic_avglamscl, rn_sic_avgphiscl, & 294 & nn_1dint, nn_2dint_default, & 295 & nn_2dint_sla, nn_2dint_sst, & 296 & nn_2dint_sss, nn_2dint_sic, & 297 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 298 & nn_profdavtypes 205 299 206 300 !----------------------------------------------------------------------- … … 208 302 !----------------------------------------------------------------------- 209 303 210 enactfiles(:) = '' 211 coriofiles(:) = '' 212 profbfiles(:) = '' 213 slafilesact(:) = '' 214 slafilespas(:) = '' 215 slafbfiles(:) = '' 216 sstfiles(:) = '' 217 sstfbfiles(:) = '' 218 seaicefiles(:) = '' 219 velcurfiles(:) = '' 220 veladcpfiles(:) = '' 221 velavcurfiles(:) = '' 222 velhrcurfiles(:) = '' 223 velavadcpfiles(:) = '' 224 velhradcpfiles(:) = '' 225 velfbfiles(:) = '' 226 velcurfiles(:) = '' 227 veladcpfiles(:) = '' 228 endailyavtypes(:) = -1 229 endailyavtypes(1) = 820 230 ln_profb_ena(:) = .FALSE. 231 ln_profb_enatim(:) = .TRUE. 232 ln_velfb_av(:) = .FALSE. 233 ln_ignmis = .FALSE. 234 235 CALL ini_date( dobsini ) 236 CALL fin_date( dobsend ) 237 238 ! Read Namelist namobs : control observation diagnostics 239 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation 304 ! Some namelist arrays need initialising 305 cn_profbfiles(:) = '' 306 cn_slafbfiles(:) = '' 307 cn_sstfbfiles(:) = '' 308 cn_sicfbfiles(:) = '' 309 cn_velfbfiles(:) = '' 310 cn_sssfbfiles(:) = '' 311 cn_slchltotfbfiles(:) = '' 312 cn_slchldiafbfiles(:) = '' 313 cn_slchlnonfbfiles(:) = '' 314 cn_slchldinfbfiles(:) = '' 315 cn_slchlmicfbfiles(:) = '' 316 cn_slchlnanfbfiles(:) = '' 317 cn_slchlpicfbfiles(:) = '' 318 cn_schltotfbfiles(:) = '' 319 cn_slphytotfbfiles(:) = '' 320 cn_slphydiafbfiles(:) = '' 321 cn_slphynonfbfiles(:) = '' 322 cn_sspmfbfiles(:) = '' 323 cn_skd490fbfiles(:) = '' 324 cn_sfco2fbfiles(:) = '' 325 cn_spco2fbfiles(:) = '' 326 cn_plchltotfbfiles(:) = '' 327 cn_pchltotfbfiles(:) = '' 328 cn_pno3fbfiles(:) = '' 329 cn_psi4fbfiles(:) = '' 330 cn_ppo4fbfiles(:) = '' 331 cn_pdicfbfiles(:) = '' 332 cn_palkfbfiles(:) = '' 333 cn_pphfbfiles(:) = '' 334 cn_po2fbfiles(:) = '' 335 cn_sstbiasfiles(:) = '' 336 nn_profdavtypes(:) = -1 337 338 CALL ini_date( rn_dobsini ) 339 CALL fin_date( rn_dobsend ) 340 341 ! Read namelist namobs : control observation diagnostics 342 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 240 343 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 241 344 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 242 345 243 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation346 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 244 347 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 245 348 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 246 349 IF(lwm) WRITE ( numond, namobs ) 247 350 248 ! Count number of files for each type 249 IF (ln_ena) THEN 250 lmask(:) = .FALSE. 251 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 252 jnumenact = COUNT(lmask) 351 lk_diaobs = .FALSE. 352 #if defined key_diaobs 353 IF ( ln_diaobs ) lk_diaobs = .TRUE. 354 #endif 355 356 IF ( .NOT. lk_diaobs ) THEN 357 IF(lwp) WRITE(numout,cform_war) 358 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' 359 RETURN 253 360 ENDIF 254 IF (ln_cor) THEN 255 lmask(:) = .FALSE. 256 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 257 jnumcorio = COUNT(lmask) 258 ENDIF 259 IF (ln_profb) THEN 260 lmask(:) = .FALSE. 261 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 262 jnumprofb = COUNT(lmask) 263 ENDIF 264 IF (ln_sladt) THEN 265 lmask(:) = .FALSE. 266 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 267 jnumslaact = COUNT(lmask) 268 lmask(:) = .FALSE. 269 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 270 jnumslapas = COUNT(lmask) 271 ENDIF 272 IF (ln_slafb) THEN 273 lmask(:) = .FALSE. 274 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 275 jnumslafb = COUNT(lmask) 276 lmask(:) = .FALSE. 277 ENDIF 278 IF (ln_ghrsst) THEN 279 lmask(:) = .FALSE. 280 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 281 jnumsst = COUNT(lmask) 282 ENDIF 283 IF (ln_sstfb) THEN 284 lmask(:) = .FALSE. 285 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 286 jnumsstfb = COUNT(lmask) 287 lmask(:) = .FALSE. 288 ENDIF 289 IF (ln_seaice) THEN 290 lmask(:) = .FALSE. 291 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 292 jnumseaice = COUNT(lmask) 293 ENDIF 294 IF (ln_velavcur) THEN 295 lmask(:) = .FALSE. 296 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 297 jnumvelavcur = COUNT(lmask) 298 ENDIF 299 IF (ln_velhrcur) THEN 300 lmask(:) = .FALSE. 301 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 302 jnumvelhrcur = COUNT(lmask) 303 ENDIF 304 IF (ln_velavadcp) THEN 305 lmask(:) = .FALSE. 306 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 307 jnumvelavadcp = COUNT(lmask) 308 ENDIF 309 IF (ln_velhradcp) THEN 310 lmask(:) = .FALSE. 311 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 312 jnumvelhradcp = COUNT(lmask) 313 ENDIF 314 IF (ln_velfb) THEN 315 lmask(:) = .FALSE. 316 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 317 jnumvelfb = COUNT(lmask) 318 lmask(:) = .FALSE. 319 ENDIF 320 321 ! Control print 361 322 362 IF(lwp) THEN 323 363 WRITE(numout,*) … … 325 365 WRITE(numout,*) '~~~~~~~~~~~~' 326 366 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 327 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 328 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 329 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena 330 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor 331 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb 332 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 333 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt 334 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb 335 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 336 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 337 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst 338 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst 339 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb 340 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 341 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 342 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 343 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 344 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 345 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 346 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 347 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 348 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 349 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 350 WRITE(numout,*) & 351 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 367 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 368 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 369 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 370 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 371 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 372 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 373 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 374 WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot 375 WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia 376 WRITE(numout,*) ' Logical switch for surface non-diatom logchl obs ln_slchlnon = ', ln_slchlnon 377 WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin 378 WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic 379 WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan 380 WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic 381 WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot 382 WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot 383 WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia 384 WRITE(numout,*) ' Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 385 WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm 386 WRITE(numout,*) ' Logical switch for surface Kd490 observations ln_skd490 = ', ln_skd490 387 WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 388 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 389 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot 390 WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot 391 WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 392 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 393 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 394 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic 395 WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk 396 WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph 397 WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 398 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 399 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 352 400 IF (ln_grid_search_lookup) & 353 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 354 IF (ln_ena) THEN 355 DO ji = 1, jnumenact 356 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 357 TRIM(enactfiles(ji)) 358 END DO 359 ENDIF 360 IF (ln_cor) THEN 361 DO ji = 1, jnumcorio 362 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 363 TRIM(coriofiles(ji)) 364 END DO 365 ENDIF 366 IF (ln_profb) THEN 367 DO ji = 1, jnumprofb 368 IF (ln_profb_ena(ji)) THEN 369 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 370 TRIM(profbfiles(ji)) 401 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 402 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 403 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 404 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 405 WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default 406 WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla 407 WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst 408 WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss 409 WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic 410 WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl 411 WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl 412 WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs 413 WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl 414 WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl 415 WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs 416 WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl 417 WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl 418 WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs 419 WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl 420 WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl 421 WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs 422 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 423 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 424 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 425 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 426 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 427 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 428 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 429 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 430 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 431 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 432 WRITE(numout,*) ' Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim 433 WRITE(numout,*) ' Logical switch for time-mean of SLA ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg 434 ENDIF 435 !----------------------------------------------------------------------- 436 ! Set up list of observation types to be used 437 ! and the files associated with each type 438 !----------------------------------------------------------------------- 439 440 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & 441 & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & 442 & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) 443 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 444 & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 445 & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & 446 & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & 447 & ln_skd490, ln_sfco2, ln_spco2 /) ) 448 449 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 450 IF(lwp) WRITE(numout,cform_war) 451 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 452 & ' are set to .FALSE. so turning off calls to dia_obs' 453 nwarn = nwarn + 1 454 lk_diaobs = .FALSE. 455 RETURN 456 ENDIF 457 458 IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN 459 IF(lwp) WRITE(numout,cform_war) 460 IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', & 461 & ' so climatological T/S not available and will not be output' 462 nwarn = nwarn + 1 463 ln_output_clim = .FALSE. 464 ENDIF 465 466 467 IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes 468 IF ( nproftypes > 0 ) THEN 469 470 ALLOCATE( cobstypesprof(nproftypes) ) 471 ALLOCATE( ifilesprof(nproftypes) ) 472 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 473 474 jtype = 0 475 IF (ln_t3d .OR. ln_s3d) THEN 476 jtype = jtype + 1 477 cobstypesprof(jtype) = 'prof' 478 clproffiles(jtype,:) = cn_profbfiles 479 ENDIF 480 IF (ln_vel3d) THEN 481 jtype = jtype + 1 482 cobstypesprof(jtype) = 'vel' 483 clproffiles(jtype,:) = cn_velfbfiles 484 ENDIF 485 IF (ln_plchltot) THEN 486 jtype = jtype + 1 487 cobstypesprof(jtype) = 'plchltot' 488 clproffiles(jtype,:) = cn_plchltotfbfiles 489 ENDIF 490 IF (ln_pchltot) THEN 491 jtype = jtype + 1 492 cobstypesprof(jtype) = 'pchltot' 493 clproffiles(jtype,:) = cn_pchltotfbfiles 494 ENDIF 495 IF (ln_pno3) THEN 496 jtype = jtype + 1 497 cobstypesprof(jtype) = 'pno3' 498 clproffiles(jtype,:) = cn_pno3fbfiles 499 ENDIF 500 IF (ln_psi4) THEN 501 jtype = jtype + 1 502 cobstypesprof(jtype) = 'psi4' 503 clproffiles(jtype,:) = cn_psi4fbfiles 504 ENDIF 505 IF (ln_ppo4) THEN 506 jtype = jtype + 1 507 cobstypesprof(jtype) = 'ppo4' 508 clproffiles(jtype,:) = cn_ppo4fbfiles 509 ENDIF 510 IF (ln_pdic) THEN 511 jtype = jtype + 1 512 cobstypesprof(jtype) = 'pdic' 513 clproffiles(jtype,:) = cn_pdicfbfiles 514 ENDIF 515 IF (ln_palk) THEN 516 jtype = jtype + 1 517 cobstypesprof(jtype) = 'palk' 518 clproffiles(jtype,:) = cn_palkfbfiles 519 ENDIF 520 IF (ln_pph) THEN 521 jtype = jtype + 1 522 cobstypesprof(jtype) = 'pph' 523 clproffiles(jtype,:) = cn_pphfbfiles 524 ENDIF 525 IF (ln_po2) THEN 526 jtype = jtype + 1 527 cobstypesprof(jtype) = 'po2' 528 clproffiles(jtype,:) = cn_po2fbfiles 529 ENDIF 530 531 CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 532 533 ENDIF 534 535 IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes 536 IF ( nsurftypes > 0 ) THEN 537 538 ALLOCATE( cobstypessurf(nsurftypes) ) 539 ALLOCATE( ifilessurf(nsurftypes) ) 540 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 541 ALLOCATE(n2dintsurf(nsurftypes)) 542 ALLOCATE(ravglamscl(nsurftypes)) 543 ALLOCATE(ravgphiscl(nsurftypes)) 544 ALLOCATE(lfpindegs(nsurftypes)) 545 ALLOCATE(llnightav(nsurftypes)) 546 547 jtype = 0 548 IF (ln_sla) THEN 549 jtype = jtype + 1 550 cobstypessurf(jtype) = 'sla' 551 clsurffiles(jtype,:) = cn_slafbfiles 552 ENDIF 553 IF (ln_sst) THEN 554 jtype = jtype + 1 555 cobstypessurf(jtype) = 'sst' 556 clsurffiles(jtype,:) = cn_sstfbfiles 557 ENDIF 558 IF (ln_sic) THEN 559 jtype = jtype + 1 560 cobstypessurf(jtype) = 'sic' 561 clsurffiles(jtype,:) = cn_sicfbfiles 562 ENDIF 563 IF (ln_sss) THEN 564 jtype = jtype + 1 565 cobstypessurf(jtype) = 'sss' 566 clsurffiles(jtype,:) = cn_sssfbfiles 567 ENDIF 568 IF (ln_slchltot) THEN 569 jtype = jtype + 1 570 cobstypessurf(jtype) = 'slchltot' 571 clsurffiles(jtype,:) = cn_slchltotfbfiles 572 ENDIF 573 IF (ln_slchldia) THEN 574 jtype = jtype + 1 575 cobstypessurf(jtype) = 'slchldia' 576 clsurffiles(jtype,:) = cn_slchldiafbfiles 577 ENDIF 578 IF (ln_slchlnon) THEN 579 jtype = jtype + 1 580 cobstypessurf(jtype) = 'slchlnon' 581 clsurffiles(jtype,:) = cn_slchlnonfbfiles 582 ENDIF 583 IF (ln_slchldin) THEN 584 jtype = jtype + 1 585 cobstypessurf(jtype) = 'slchldin' 586 clsurffiles(jtype,:) = cn_slchldinfbfiles 587 ENDIF 588 IF (ln_slchlmic) THEN 589 jtype = jtype + 1 590 cobstypessurf(jtype) = 'slchlmic' 591 clsurffiles(jtype,:) = cn_slchlmicfbfiles 592 ENDIF 593 IF (ln_slchlnan) THEN 594 jtype = jtype + 1 595 cobstypessurf(jtype) = 'slchlnan' 596 clsurffiles(jtype,:) = cn_slchlnanfbfiles 597 ENDIF 598 IF (ln_slchlpic) THEN 599 jtype = jtype + 1 600 cobstypessurf(jtype) = 'slchlpic' 601 clsurffiles(jtype,:) = cn_slchlpicfbfiles 602 ENDIF 603 IF (ln_schltot) THEN 604 jtype = jtype + 1 605 cobstypessurf(jtype) = 'schltot' 606 clsurffiles(jtype,:) = cn_schltotfbfiles 607 ENDIF 608 IF (ln_slphytot) THEN 609 jtype = jtype + 1 610 cobstypessurf(jtype) = 'slphytot' 611 clsurffiles(jtype,:) = cn_slphytotfbfiles 612 ENDIF 613 IF (ln_slphydia) THEN 614 jtype = jtype + 1 615 cobstypessurf(jtype) = 'slphydia' 616 clsurffiles(jtype,:) = cn_slphydiafbfiles 617 ENDIF 618 IF (ln_slphynon) THEN 619 jtype = jtype + 1 620 cobstypessurf(jtype) = 'slphynon' 621 clsurffiles(jtype,:) = cn_slphynonfbfiles 622 ENDIF 623 IF (ln_sspm) THEN 624 jtype = jtype + 1 625 cobstypessurf(jtype) = 'sspm' 626 clsurffiles(jtype,:) = cn_sspmfbfiles 627 ENDIF 628 IF (ln_skd490) THEN 629 jtype = jtype + 1 630 cobstypessurf(jtype) = 'skd490' 631 clsurffiles(jtype,:) = cn_skd490fbfiles 632 ENDIF 633 IF (ln_sfco2) THEN 634 jtype = jtype + 1 635 cobstypessurf(jtype) = 'sfco2' 636 clsurffiles(jtype,:) = cn_sfco2fbfiles 637 ENDIF 638 IF (ln_spco2) THEN 639 jtype = jtype + 1 640 cobstypessurf(jtype) = 'spco2' 641 clsurffiles(jtype,:) = cn_spco2fbfiles 642 ENDIF 643 644 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 645 646 DO jtype = 1, nsurftypes 647 648 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 649 IF ( nn_2dint_sla == -1 ) THEN 650 n2dint_type = nn_2dint_default 371 651 ELSE 372 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 373 TRIM(profbfiles(ji)) 652 n2dint_type = nn_2dint_sla 374 653 ENDIF 375 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 376 END DO 377 ENDIF 378 IF (ln_sladt) THEN 379 DO ji = 1, jnumslaact 380 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 381 TRIM(slafilesact(ji)) 382 END DO 383 DO ji = 1, jnumslapas 384 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 385 TRIM(slafilespas(ji)) 386 END DO 387 ENDIF 388 IF (ln_slafb) THEN 389 DO ji = 1, jnumslafb 390 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 391 TRIM(slafbfiles(ji)) 392 END DO 393 ENDIF 394 IF (ln_ghrsst) THEN 395 DO ji = 1, jnumsst 396 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 397 TRIM(sstfiles(ji)) 398 END DO 399 ENDIF 400 IF (ln_sstfb) THEN 401 DO ji = 1, jnumsstfb 402 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 403 TRIM(sstfbfiles(ji)) 404 END DO 405 ENDIF 406 IF (ln_seaice) THEN 407 DO ji = 1, jnumseaice 408 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 409 TRIM(seaicefiles(ji)) 410 END DO 411 ENDIF 412 IF (ln_velavcur) THEN 413 DO ji = 1, jnumvelavcur 414 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 415 TRIM(velavcurfiles(ji)) 416 END DO 417 ENDIF 418 IF (ln_velhrcur) THEN 419 DO ji = 1, jnumvelhrcur 420 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 421 TRIM(velhrcurfiles(ji)) 422 END DO 423 ENDIF 424 IF (ln_velavadcp) THEN 425 DO ji = 1, jnumvelavadcp 426 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 427 TRIM(velavadcpfiles(ji)) 428 END DO 429 ENDIF 430 IF (ln_velhradcp) THEN 431 DO ji = 1, jnumvelhradcp 432 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 433 TRIM(velhradcpfiles(ji)) 434 END DO 435 ENDIF 436 IF (ln_velfb) THEN 437 DO ji = 1, jnumvelfb 438 IF (ln_velfb_av(ji)) THEN 439 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 440 TRIM(velfbfiles(ji)) 654 ztype_avglamscl = rn_sla_avglamscl 655 ztype_avgphiscl = rn_sla_avgphiscl 656 ltype_fp_indegs = ln_sla_fp_indegs 657 ltype_night = .FALSE. 658 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 659 IF ( nn_2dint_sst == -1 ) THEN 660 n2dint_type = nn_2dint_default 441 661 ELSE 442 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 443 TRIM(velfbfiles(ji)) 662 n2dint_type = nn_2dint_sst 444 663 ENDIF 445 END DO 446 ENDIF 447 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 448 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 449 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint 450 WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint 451 WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea 452 WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc 453 WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr 454 WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff 455 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 456 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 457 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 664 ztype_avglamscl = rn_sst_avglamscl 665 ztype_avgphiscl = rn_sst_avgphiscl 666 ltype_fp_indegs = ln_sst_fp_indegs 667 ltype_night = ln_sstnight 668 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 669 IF ( nn_2dint_sic == -1 ) THEN 670 n2dint_type = nn_2dint_default 671 ELSE 672 n2dint_type = nn_2dint_sic 673 ENDIF 674 ztype_avglamscl = rn_sic_avglamscl 675 ztype_avgphiscl = rn_sic_avgphiscl 676 ltype_fp_indegs = ln_sic_fp_indegs 677 ltype_night = .FALSE. 678 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 679 IF ( nn_2dint_sss == -1 ) THEN 680 n2dint_type = nn_2dint_default 681 ELSE 682 n2dint_type = nn_2dint_sss 683 ENDIF 684 ztype_avglamscl = rn_sss_avglamscl 685 ztype_avgphiscl = rn_sss_avgphiscl 686 ltype_fp_indegs = ln_sss_fp_indegs 687 ltype_night = .FALSE. 688 ELSE 689 n2dint_type = nn_2dint_default 690 ztype_avglamscl = rn_default_avglamscl 691 ztype_avgphiscl = rn_default_avgphiscl 692 ltype_fp_indegs = ln_default_fp_indegs 693 ltype_night = .FALSE. 694 ENDIF 695 696 CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 697 & nn_2dint_default, n2dint_type, & 698 & ztype_avglamscl, ztype_avgphiscl, & 699 & ltype_fp_indegs, ltype_night, & 700 & n2dintsurf, ravglamscl, ravgphiscl, & 701 & lfpindegs, llnightav ) 702 703 END DO 458 704 459 705 ENDIF 460 706 707 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 708 709 710 !----------------------------------------------------------------------- 711 ! Obs operator parameter checking and initialisations 712 !----------------------------------------------------------------------- 713 461 714 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 462 715 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 464 717 ENDIF 465 718 466 CALL obs_typ_init 467 468 CALL mppmap_init 469 470 ! Parameter control 471 #if defined key_diaobs 472 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 473 & ( .NOT. ln_vel3d ).AND. & 474 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 475 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 476 IF(lwp) WRITE(numout,cform_war) 477 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 478 & ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 479 nwarn = nwarn + 1 480 ENDIF 481 #endif 482 483 CALL obs_grid_setup( ) 484 IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 719 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 485 720 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 486 721 & ' is not available') 487 722 ENDIF 488 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 489 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 723 724 IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN 725 CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & 490 726 & ' is not available') 491 727 ENDIF 728 729 CALL obs_typ_init 730 731 CALL mppmap_init 732 733 CALL obs_grid_setup( ) 492 734 493 735 !----------------------------------------------------------------------- 494 736 ! Depending on switches read the various observation types 495 737 !----------------------------------------------------------------------- 496 ! - Temperature/salinity profiles 497 498 IF ( ln_t3d .OR. ln_s3d ) THEN 499 500 ! Set the number of variables for profiles to 2 (T and S) 501 nprofvars = 2 502 ! Set the number of extra variables for profiles to 1 (insitu temp). 503 nprofextr = 1 504 505 ! Count how may insitu data sets we have and allocate data. 506 jprofset = 0 507 IF ( ln_ena ) jprofset = jprofset + 1 508 IF ( ln_cor ) jprofset = jprofset + 1 509 IF ( ln_profb ) jprofset = jprofset + jnumprofb 510 nprofsets = jprofset 511 IF ( nprofsets > 0 ) THEN 512 ALLOCATE(ld_enact(nprofsets)) 513 ALLOCATE(profdata(nprofsets)) 514 ALLOCATE(prodatqc(nprofsets)) 515 ENDIF 516 517 jprofset = 0 518 519 ! ENACT insitu data 520 521 IF ( ln_ena ) THEN 522 523 jprofset = jprofset + 1 738 739 IF ( nproftypes > 0 ) THEN 740 741 ALLOCATE(profdata(nproftypes)) 742 ALLOCATE(profdataqc(nproftypes)) 743 ALLOCATE(nvarsprof(nproftypes)) 744 ALLOCATE(nextrprof(nproftypes)) 745 746 DO jtype = 1, nproftypes 524 747 525 ld_enact(jprofset) = .TRUE. 526 527 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 528 & jnumenact, enactfiles(1:jnumenact), & 529 & nprofvars, nprofextr, & 530 & nitend-nit000+2, & 531 & dobsini, dobsend, ln_t3d, ln_s3d, & 532 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 533 & kdailyavtypes = endailyavtypes ) 534 535 DO jvar = 1, 2 536 537 CALL obs_prof_staend( profdata(jprofset), jvar ) 538 748 ltype_clim = .FALSE. 749 750 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 751 nvarsprof(jtype) = 2 752 nextrprof(jtype) = 1 753 IF ( ln_output_clim ) ltype_clim = .TRUE. 754 ALLOCATE(llvar(nvarsprof(jtype))) 755 ALLOCATE(clvars(nvarsprof(jtype))) 756 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 757 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 758 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 759 llvar(1) = ln_t3d 760 llvar(2) = ln_s3d 761 clvars(1) = 'POTM' 762 clvars(2) = 'PSAL' 763 zglam(:,:,1) = glamt(:,:) 764 zglam(:,:,2) = glamt(:,:) 765 zgphi(:,:,1) = gphit(:,:) 766 zgphi(:,:,2) = gphit(:,:) 767 zmask(:,:,:,1) = tmask(:,:,:) 768 zmask(:,:,:,2) = tmask(:,:,:) 769 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 770 nvarsprof(jtype) = 2 771 nextrprof(jtype) = 2 772 ALLOCATE(llvar(nvarsprof(jtype))) 773 ALLOCATE(clvars(nvarsprof(jtype))) 774 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 775 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 776 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 777 llvar(1) = ln_vel3d 778 llvar(2) = ln_vel3d 779 clvars(1) = 'UVEL' 780 clvars(2) = 'VVEL' 781 zglam(:,:,1) = glamu(:,:) 782 zglam(:,:,2) = glamv(:,:) 783 zgphi(:,:,1) = gphiu(:,:) 784 zgphi(:,:,2) = gphiv(:,:) 785 zmask(:,:,:,1) = umask(:,:,:) 786 zmask(:,:,:,2) = vmask(:,:,:) 787 ELSE 788 nvarsprof(jtype) = 1 789 nextrprof(jtype) = 0 790 ALLOCATE(llvar(nvarsprof(jtype))) 791 ALLOCATE(clvars(nvarsprof(jtype))) 792 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 793 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 794 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 795 llvar(1) = .TRUE. 796 zglam(:,:,1) = glamt(:,:) 797 zgphi(:,:,1) = gphit(:,:) 798 zmask(:,:,:,1) = tmask(:,:,:) 799 IF ( TRIM(cobstypesprof(jtype)) == 'plchltot' ) THEN 800 clvars(1) = 'PLCHLTOT' 801 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pchltot' ) THEN 802 clvars(1) = 'PCHLTOT' 803 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pno3' ) THEN 804 clvars(1) = 'PNO3' 805 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'psi4' ) THEN 806 clvars(1) = 'PSI4' 807 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'ppo4' ) THEN 808 clvars(1) = 'PPO4' 809 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pdic' ) THEN 810 clvars(1) = 'PDIC' 811 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'palk' ) THEN 812 clvars(1) = 'PALK' 813 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pph' ) THEN 814 clvars(1) = 'PPH' 815 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'po2' ) THEN 816 clvars(1) = 'PO2' 817 ENDIF 818 ENDIF 819 820 !Read in profile or profile obs types 821 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 822 & clproffiles(jtype,1:ifilesprof(jtype)), & 823 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 824 & rn_dobsini, rn_dobsend, llvar, & 825 & ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, clvars, & 826 & kdailyavtypes = nn_profdavtypes ) 827 828 DO jvar = 1, nvarsprof(jtype) 829 CALL obs_prof_staend( profdata(jtype), jvar ) 539 830 END DO 540 831 541 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 542 & ln_t3d, ln_s3d, ln_nea, & 543 & kdailyavtypes=endailyavtypes ) 832 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 833 & llvar, & 834 & jpi, jpj, jpk, & 835 & zmask, zglam, zgphi, & 836 & ln_nea, ln_bound_reject, & 837 & kdailyavtypes = nn_profdavtypes ) 544 838 545 ENDIF 546 547 ! Coriolis insitu data 548 549 IF ( ln_cor ) THEN 550 551 jprofset = jprofset + 1 552 553 ld_enact(jprofset) = .FALSE. 554 555 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 556 & jnumcorio, coriofiles(1:jnumcorio), & 557 & nprofvars, nprofextr, & 558 & nitend-nit000+2, & 559 & dobsini, dobsend, ln_t3d, ln_s3d, & 560 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 561 562 DO jvar = 1, 2 563 564 CALL obs_prof_staend( profdata(jprofset), jvar ) 565 566 END DO 567 568 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 569 & ln_t3d, ln_s3d, ln_nea ) 839 DEALLOCATE( llvar, clvars ) 840 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zglam ) 841 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zgphi ) 842 CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 843 844 END DO 845 846 DEALLOCATE( ifilesprof, clproffiles ) 847 848 ENDIF 849 850 IF ( nsurftypes > 0 ) THEN 851 852 ALLOCATE(surfdata(nsurftypes)) 853 ALLOCATE(surfdataqc(nsurftypes)) 854 ALLOCATE(nvarssurf(nsurftypes)) 855 ALLOCATE(nextrsurf(nsurftypes)) 856 857 DO jtype = 1, nsurftypes 858 859 ltype_clim = .FALSE. 860 IF ( ln_output_clim .AND. & 861 & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. & 862 & ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) & 863 & ltype_clim = .TRUE. 864 865 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 866 nvarssurf(jtype) = 1 867 nextrsurf(jtype) = 2 868 ELSE 869 nvarssurf(jtype) = 1 870 nextrsurf(jtype) = 0 871 ENDIF 570 872 571 ENDIF 572 573 ! Feedback insitu data 574 575 IF ( ln_profb ) THEN 576 577 DO jset = 1, jnumprofb 578 579 jprofset = jprofset + 1 580 ld_enact (jprofset) = ln_profb_ena(jset) 581 582 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 583 & 1, profbfiles(jset:jset), & 584 & nprofvars, nprofextr, & 585 & nitend-nit000+2, & 586 & dobsini, dobsend, ln_t3d, ln_s3d, & 587 & ln_ignmis, ln_s_at_t, & 588 & ld_enact(jprofset).AND.& 589 & ln_profb_enatim(jset), & 590 & .FALSE., kdailyavtypes = endailyavtypes ) 591 592 DO jvar = 1, 2 593 594 CALL obs_prof_staend( profdata(jprofset), jvar ) 595 873 ALLOCATE( clvars( nvarssurf(jtype) ) ) 874 875 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 876 clvars(1) = 'SLA' 877 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 878 clvars(1) = 'SST' 879 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 880 clvars(1) = 'ICECONC' 881 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 882 clvars(1) = 'SSS' 883 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN 884 clvars(1) = 'SLCHLTOT' 885 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldia' ) THEN 886 clvars(1) = 'SLCHLDIA' 887 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnon' ) THEN 888 clvars(1) = 'SLCHLNON' 889 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldin' ) THEN 890 clvars(1) = 'SLCHLDIN' 891 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlmic' ) THEN 892 clvars(1) = 'SLCHLMIC' 893 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnan' ) THEN 894 clvars(1) = 'SLCHLNAN' 895 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlpic' ) THEN 896 clvars(1) = 'SLCHLPIC' 897 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'schltot' ) THEN 898 clvars(1) = 'SCHLTOT' 899 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphytot' ) THEN 900 clvars(1) = 'SLPHYTOT' 901 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphydia' ) THEN 902 clvars(1) = 'SLPHYDIA' 903 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphynon' ) THEN 904 clvars(1) = 'SLPHYNON' 905 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sspm' ) THEN 906 clvars(1) = 'SSPM' 907 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'skd490' ) THEN 908 clvars(1) = 'SKD490' 909 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sfco2' ) THEN 910 clvars(1) = 'SFCO2' 911 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'spco2' ) THEN 912 clvars(1) = 'SPCO2' 913 ENDIF 914 915 !Read in surface obs types 916 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 917 & clsurffiles(jtype,1:ifilessurf(jtype)), & 918 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 919 & rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., & 920 & llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 921 922 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 923 924 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 925 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 926 IF ( ln_altbias ) & 927 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 928 ENDIF 929 930 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 931 jnumsstbias = 0 932 DO jfile = 1, jpmaxnfiles 933 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 934 & jnumsstbias = jnumsstbias + 1 596 935 END DO 597 598 IF ( ld_enact(jprofset) ) THEN 599 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 600 & ln_t3d, ln_s3d, ln_nea, & 601 & kdailyavtypes = endailyavtypes ) 602 ELSE 603 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 604 & ln_t3d, ln_s3d, ln_nea ) 936 IF ( jnumsstbias == 0 ) THEN 937 CALL ctl_stop("ln_sstbias set but no bias files to read in") 605 938 ENDIF 606 607 END DO 608 609 ENDIF 939 940 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 941 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 942 943 ENDIF 944 945 DEALLOCATE( clvars ) 946 947 END DO 948 949 DEALLOCATE( ifilessurf, clsurffiles ) 610 950 611 951 ENDIF 612 952 613 ! - Sea level anomalies614 IF ( ln_sla ) THEN615 ! Set the number of variables for sla to 1616 nslavars = 1617 618 ! Set the number of extra variables for sla to 2619 nslaextr = 2620 621 ! Set the number of sla data sets to 2622 nslasets = 0623 IF ( ln_sladt ) THEN624 nslasets = nslasets + 2625 ENDIF626 IF ( ln_slafb ) THEN627 nslasets = nslasets + jnumslafb628 ENDIF629 630 ALLOCATE(sladata(nslasets))631 ALLOCATE(sladatqc(nslasets))632 sladata(:)%nsurf=0633 sladatqc(:)%nsurf=0634 635 nslasets = 0636 637 ! AVISO SLA data638 639 IF ( ln_sladt ) THEN640 641 ! Active SLA observations642 643 nslasets = nslasets + 1644 645 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &646 & slafilesact(1:jnumslaact), &647 & nslavars, nslaextr, nitend-nit000+2, &648 & dobsini, dobsend, ln_ignmis, .FALSE. )649 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &650 & ln_sla, ln_nea )651 652 ! Passive SLA observations653 654 nslasets = nslasets + 1655 656 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &657 & slafilespas(1:jnumslapas), &658 & nslavars, nslaextr, nitend-nit000+2, &659 & dobsini, dobsend, ln_ignmis, .FALSE. )660 661 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &662 & ln_sla, ln_nea )663 664 ENDIF665 666 ! Feedback SLA data667 668 IF ( ln_slafb ) THEN669 670 DO jset = 1, jnumslafb671 672 nslasets = nslasets + 1673 674 CALL obs_rea_sla( 0, sladata(nslasets), 1, &675 & slafbfiles(jset:jset), &676 & nslavars, nslaextr, nitend-nit000+2, &677 & dobsini, dobsend, ln_ignmis, .FALSE. )678 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &679 & ln_sla, ln_nea )680 681 END DO682 683 ENDIF684 685 CALL obs_rea_mdt( nslasets, sladatqc, n2dint )686 687 ! read in altimeter bias688 689 IF ( ln_altbias ) THEN690 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )691 ENDIF692 693 ENDIF694 695 ! - Sea surface height696 IF ( ln_ssh ) THEN697 IF(lwp) WRITE(numout,*) ' SSH currently not available'698 ENDIF699 700 ! - Sea surface temperature701 IF ( ln_sst ) THEN702 703 ! Set the number of variables for sst to 1704 nsstvars = 1705 706 ! Set the number of extra variables for sst to 0707 nsstextr = 0708 709 nsstsets = 0710 711 IF (ln_reysst) nsstsets = nsstsets + 1712 IF (ln_ghrsst) nsstsets = nsstsets + 1713 IF ( ln_sstfb ) THEN714 nsstsets = nsstsets + jnumsstfb715 ENDIF716 717 ALLOCATE(sstdata(nsstsets))718 ALLOCATE(sstdatqc(nsstsets))719 ALLOCATE(ld_sstnight(nsstsets))720 sstdata(:)%nsurf=0721 sstdatqc(:)%nsurf=0722 ld_sstnight(:)=.false.723 724 nsstsets = 0725 726 IF (ln_reysst) THEN727 728 nsstsets = nsstsets + 1729 730 ld_sstnight(nsstsets) = ln_sstnight731 732 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &733 & nsstvars, nsstextr, &734 & nitend-nit000+2, dobsini, dobsend )735 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &736 & ln_nea )737 738 ENDIF739 740 IF (ln_ghrsst) THEN741 742 nsstsets = nsstsets + 1743 744 ld_sstnight(nsstsets) = ln_sstnight745 746 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &747 & sstfiles(1:jnumsst), &748 & nsstvars, nsstextr, nitend-nit000+2, &749 & dobsini, dobsend, ln_ignmis, .FALSE. )750 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &751 & ln_nea )752 753 ENDIF754 755 ! Feedback SST data756 757 IF ( ln_sstfb ) THEN758 759 DO jset = 1, jnumsstfb760 761 nsstsets = nsstsets + 1762 763 ld_sstnight(nsstsets) = ln_sstnight764 765 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &766 & sstfbfiles(jset:jset), &767 & nsstvars, nsstextr, nitend-nit000+2, &768 & dobsini, dobsend, ln_ignmis, .FALSE. )769 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &770 & ln_sst, ln_nea )771 772 END DO773 774 ENDIF775 776 ENDIF777 778 ! - Sea surface salinity779 IF ( ln_sss ) THEN780 IF(lwp) WRITE(numout,*) ' SSS currently not available'781 ENDIF782 783 ! - Sea Ice Concentration784 785 IF ( ln_seaice ) THEN786 787 ! Set the number of variables for seaice to 1788 nseaicevars = 1789 790 ! Set the number of extra variables for seaice to 0791 nseaiceextr = 0792 793 ! Set the number of data sets to 1794 nseaicesets = 1795 796 ALLOCATE(seaicedata(nseaicesets))797 ALLOCATE(seaicedatqc(nseaicesets))798 seaicedata(:)%nsurf=0799 seaicedatqc(:)%nsurf=0800 801 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &802 & seaicefiles(1:jnumseaice), &803 & nseaicevars, nseaiceextr, nitend-nit000+2, &804 & dobsini, dobsend, ln_ignmis, .FALSE. )805 806 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &807 & ln_seaice, ln_nea )808 809 ENDIF810 811 IF (ln_vel3d) THEN812 813 ! Set the number of variables for profiles to 2 (U and V)814 nvelovars = 2815 816 ! Set the number of extra variables for profiles to 2 to store817 ! rotation parameters818 nveloextr = 2819 820 jveloset = 0821 822 IF ( ln_velavcur ) jveloset = jveloset + 1823 IF ( ln_velhrcur ) jveloset = jveloset + 1824 IF ( ln_velavadcp ) jveloset = jveloset + 1825 IF ( ln_velhradcp ) jveloset = jveloset + 1826 IF (ln_velfb) jveloset = jveloset + jnumvelfb827 828 nvelosets = jveloset829 IF ( nvelosets > 0 ) THEN830 ALLOCATE( velodata(nvelosets) )831 ALLOCATE( veldatqc(nvelosets) )832 ALLOCATE( ld_velav(nvelosets) )833 ENDIF834 835 jveloset = 0836 837 ! Daily averaged data838 839 IF ( ln_velavcur ) THEN840 841 jveloset = jveloset + 1842 843 ld_velav(jveloset) = .TRUE.844 845 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &846 & velavcurfiles(1:jnumvelavcur), &847 & nvelovars, nveloextr, &848 & nitend-nit000+2, &849 & dobsini, dobsend, ln_ignmis, &850 & ld_velav(jveloset), &851 & .FALSE. )852 853 DO jvar = 1, 2854 CALL obs_prof_staend( velodata(jveloset), jvar )855 END DO856 857 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &858 & ln_vel3d, ln_nea, ld_velav(jveloset) )859 860 ENDIF861 862 ! High frequency data863 864 IF ( ln_velhrcur ) THEN865 866 jveloset = jveloset + 1867 868 ld_velav(jveloset) = .FALSE.869 870 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &871 & velhrcurfiles(1:jnumvelhrcur), &872 & nvelovars, nveloextr, &873 & nitend-nit000+2, &874 & dobsini, dobsend, ln_ignmis, &875 & ld_velav(jveloset), &876 & .FALSE. )877 878 DO jvar = 1, 2879 CALL obs_prof_staend( velodata(jveloset), jvar )880 END DO881 882 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &883 & ln_vel3d, ln_nea, ld_velav(jveloset) )884 885 ENDIF886 887 ! Daily averaged data888 889 IF ( ln_velavadcp ) THEN890 891 jveloset = jveloset + 1892 893 ld_velav(jveloset) = .TRUE.894 895 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &896 & velavadcpfiles(1:jnumvelavadcp), &897 & nvelovars, nveloextr, &898 & nitend-nit000+2, &899 & dobsini, dobsend, ln_ignmis, &900 & ld_velav(jveloset), &901 & .FALSE. )902 903 DO jvar = 1, 2904 CALL obs_prof_staend( velodata(jveloset), jvar )905 END DO906 907 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &908 & ln_vel3d, ln_nea, ld_velav(jveloset) )909 910 ENDIF911 912 ! High frequency data913 914 IF ( ln_velhradcp ) THEN915 916 jveloset = jveloset + 1917 918 ld_velav(jveloset) = .FALSE.919 920 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &921 & velhradcpfiles(1:jnumvelhradcp), &922 & nvelovars, nveloextr, &923 & nitend-nit000+2, &924 & dobsini, dobsend, ln_ignmis, &925 & ld_velav(jveloset), &926 & .FALSE. )927 928 DO jvar = 1, 2929 CALL obs_prof_staend( velodata(jveloset), jvar )930 END DO931 932 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &933 & ln_vel3d, ln_nea, ld_velav(jveloset) )934 935 ENDIF936 937 IF ( ln_velfb ) THEN938 939 DO jset = 1, jnumvelfb940 941 jveloset = jveloset + 1942 943 ld_velav(jveloset) = ln_velfb_av(jset)944 945 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &946 & velfbfiles(jset:jset), &947 & nvelovars, nveloextr, &948 & nitend-nit000+2, &949 & dobsini, dobsend, ln_ignmis, &950 & ld_velav(jveloset), &951 & .FALSE. )952 953 DO jvar = 1, 2954 CALL obs_prof_staend( velodata(jveloset), jvar )955 END DO956 957 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &958 & ln_vel3d, ln_nea, ld_velav(jveloset) )959 960 961 END DO962 963 ENDIF964 965 ENDIF966 967 953 END SUBROUTINE dia_obs_init 968 954 … … 974 960 !! 975 961 !! ** Method : Call the observation operators on each time step to 976 !! compute the model equivalent of the following date: 977 !! - T profiles 978 !! - S profiles 979 !! - Sea surface height (referenced to a mean) 980 !! - Sea surface temperature 981 !! - Sea surface salinity 982 !! - Velocity component (U,V) profiles 962 !! compute the model equivalent of the following data: 963 !! - Profile data, currently T/S or U/V 964 !! - Surface data, currently SST, SLA or sea-ice concentration. 983 965 !! 984 !! ** Action : 966 !! ** Action : 985 967 !! 986 968 !! History : … … 991 973 !! ! 07-04 (G. Smith) Generalized surface operators 992 974 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 975 !! ! 15-08 (M. Martin) Combined surface/profile routines. 993 976 !!---------------------------------------------------------------------- 994 977 !! * Modules used 995 USE dom_oce, ONLY : & ! Ocean space and time domain variables996 & rdt, & 997 & gdept_1d, &998 & tmask, umask, vmask 999 USE phycst, ONLY : & ! Physical constants1000 & rday1001 USE oce, ONLY : & ! Ocean dynamics and tracers variables1002 & tsn, &1003 & un, vn,&978 USE phycst, ONLY : & ! Physical constants 979 #if defined key_fabm 980 & rt0, & 981 #endif 982 & rday 983 USE oce, ONLY : & ! Ocean dynamics and tracers variables 984 & tsn, & 985 & un, & 986 & vn, & 1004 987 & sshn 1005 988 #if defined key_lim3 1006 USE ice, ONLY : & ! LIMIce model variables989 USE ice, ONLY : & ! LIM3 Ice model variables 1007 990 & frld 1008 991 #endif 1009 992 #if defined key_lim2 1010 USE ice_2, ONLY : & ! LIMIce model variables993 USE ice_2, ONLY : & ! LIM2 Ice model variables 1011 994 & frld 1012 995 #endif 996 #if defined key_cice 997 USE sbc_oce, ONLY : fr_i ! ice fraction 998 #endif 999 #if defined key_top 1000 USE trc, ONLY : & ! Biogeochemical state variables 1001 & trn 1002 #endif 1003 #if defined key_hadocc 1004 USE par_hadocc ! HadOCC parameters 1005 USE trc, ONLY : & 1006 & HADOCC_CHL, & 1007 & HADOCC_FCO2, & 1008 & HADOCC_PCO2, & 1009 & HADOCC_FILL_FLT 1010 USE had_bgc_const, ONLY: c2n_p 1011 #elif defined key_medusa 1012 USE par_medusa ! MEDUSA parameters 1013 USE sms_medusa, ONLY: & 1014 & xthetapn, & 1015 & xthetapd 1016 #if defined key_roam 1017 USE sms_medusa, ONLY: & 1018 & f2_pco2w, & 1019 & f2_fco2w, & 1020 & f3_pH 1021 #endif 1022 #elif defined key_fabm 1023 USE par_fabm ! FABM parameters 1024 USE fabm, ONLY: & 1025 & fabm_get_interior_diagnostic_data 1026 #endif 1027 #if defined key_spm 1028 USE par_spm, ONLY: & ! Sediment parameters 1029 & jp_spm 1030 #endif 1031 1013 1032 IMPLICIT NONE 1014 1033 1015 1034 !! * Arguments 1016 INTEGER, INTENT(IN) :: kstp 1035 INTEGER, INTENT(IN) :: kstp ! Current timestep 1017 1036 !! * Local declarations 1018 INTEGER :: idaystp ! Number of timesteps per day 1019 INTEGER :: jprofset ! Profile data set loop variable 1020 INTEGER :: jslaset ! SLA data set loop variable 1021 INTEGER :: jsstset ! SST data set loop variable 1022 INTEGER :: jseaiceset ! sea ice data set loop variable 1023 INTEGER :: jveloset ! velocity profile data loop variable 1024 INTEGER :: jvar ! Variable number 1025 #if ! defined key_lim2 && ! defined key_lim3 1026 REAL(wp), POINTER, DIMENSION(:,:) :: frld 1027 #endif 1028 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1029 1030 #if ! defined key_lim2 && ! defined key_lim3 1031 CALL wrk_alloc(jpi,jpj,frld) 1032 #endif 1033 1037 INTEGER :: idaystp ! Number of timesteps per day 1038 INTEGER :: imeanstp ! Number of timesteps for sla averaging 1039 INTEGER :: jtype ! Data loop variable 1040 INTEGER :: jvar ! Variable number 1041 INTEGER :: ji, jj, jk ! Loop counters 1042 REAL(wp) :: tiny ! small number 1043 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 1044 & zprofvar, & ! Model values for variables in a prof ob 1045 & zprofclim ! Climatology values for variables in a prof ob 1046 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 1047 & zprofmask ! Mask associated with zprofvar 1048 REAL(wp), POINTER, DIMENSION(:,:) :: & 1049 & zsurfvar, & ! Model values equivalent to surface ob. 1050 & zsurfclim, & ! Climatology values for variables in a surface ob. 1051 & zsurfmask ! Mask associated with surface variable 1052 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 1053 & zglam, & ! Model longitudes for prof variables 1054 & zgphi ! Model latitudes for prof variables 1055 LOGICAL :: llog10 ! Perform log10 transform of variable 1056 #if defined key_fabm 1057 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 1058 & fabm_3d ! 3D variable from FABM 1059 #endif 1060 1034 1061 IF(lwp) THEN 1035 1062 WRITE(numout,*) 1036 1063 WRITE(numout,*) 'dia_obs : Call the observation operators', kstp 1037 1064 WRITE(numout,*) '~~~~~~~' 1065 CALL FLUSH(numout) 1038 1066 ENDIF 1039 1067 … … 1041 1069 1042 1070 !----------------------------------------------------------------------- 1043 ! No LIM => frld == 0.0_wp1071 ! Call the profile and surface observation operators 1044 1072 !----------------------------------------------------------------------- 1045 #if ! defined key_lim2 && ! defined key_lim3 1046 frld(:,:) = 0.0_wp 1047 #endif 1048 !----------------------------------------------------------------------- 1049 ! Depending on switches call various observation operators 1050 !----------------------------------------------------------------------- 1051 1052 ! - Temperature/salinity profiles 1053 IF ( ln_t3d .OR. ln_s3d ) THEN 1054 DO jprofset = 1, nprofsets 1055 IF ( ld_enact(jprofset) ) THEN 1056 CALL obs_pro_opt( prodatqc(jprofset), & 1057 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1058 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1059 & gdept_1d, tmask, n1dint, n2dint, & 1060 & kdailyavtypes = endailyavtypes ) 1073 1074 IF ( nproftypes > 0 ) THEN 1075 1076 DO jtype = 1, nproftypes 1077 1078 ! Allocate local work arrays 1079 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1080 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1081 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1082 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1083 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim ) 1084 1085 ! Defaults which might change 1086 DO jvar = 1, profdataqc(jtype)%nvar 1087 zprofmask(:,:,:,jvar) = tmask(:,:,:) 1088 zglam(:,:,jvar) = glamt(:,:) 1089 zgphi(:,:,jvar) = gphit(:,:) 1090 zprofclim(:,:,:,jvar) = 0._wp 1091 END DO 1092 1093 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 1094 1095 CASE('prof') 1096 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 1097 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 1098 IF ( ln_output_clim ) THEN 1099 zprofclim(:,:,:,1) = tclim(:,:,:) 1100 zprofclim(:,:,:,2) = sclim(:,:,:) 1101 ENDIF 1102 1103 CASE('vel') 1104 zprofvar(:,:,:,1) = un(:,:,:) 1105 zprofvar(:,:,:,2) = vn(:,:,:) 1106 zprofmask(:,:,:,1) = umask(:,:,:) 1107 zprofmask(:,:,:,2) = vmask(:,:,:) 1108 zglam(:,:,1) = glamu(:,:) 1109 zglam(:,:,2) = glamv(:,:) 1110 zgphi(:,:,1) = gphiu(:,:) 1111 zgphi(:,:,2) = gphiv(:,:) 1112 1113 CASE('plchltot') 1114 #if defined key_hadocc 1115 ! Chlorophyll from HadOCC 1116 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1117 #elif defined key_medusa 1118 ! Add non-diatom and diatom chlorophyll from MEDUSA 1119 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1120 #elif defined key_fabm 1121 ! Add all chlorophyll groups from ERSEM 1122 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1123 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1124 #else 1125 CALL ctl_stop( ' Trying to run plchltot observation operator', & 1126 & ' but no biogeochemical model appears to have been defined' ) 1127 #endif 1128 ! Take the log10 where we can, otherwise exclude 1129 tiny = 1.0e-20 1130 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 1131 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 1132 ELSEWHERE 1133 zprofvar(:,:,:,:) = obfillflt 1134 zprofmask(:,:,:,:) = 0 1135 END WHERE 1136 ! Mask out model below any excluded values, 1137 ! to avoid interpolation issues 1138 DO jvar = 1, profdataqc(jtype)%nvar 1139 DO jj = 1, jpj 1140 DO ji = 1, jpi 1141 depth_loop: DO jk = 1, jpk 1142 IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 1143 zprofmask(ji,jj,jk:jpk,jvar) = 0 1144 EXIT depth_loop 1145 ENDIF 1146 END DO depth_loop 1147 END DO 1148 END DO 1149 END DO 1150 1151 CASE('pchltot') 1152 #if defined key_hadocc 1153 ! Chlorophyll from HadOCC 1154 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1155 #elif defined key_medusa 1156 ! Add non-diatom and diatom chlorophyll from MEDUSA 1157 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1158 #elif defined key_fabm 1159 ! Add all chlorophyll groups from ERSEM 1160 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1161 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1162 #else 1163 CALL ctl_stop( ' Trying to run pchltot observation operator', & 1164 & ' but no biogeochemical model appears to have been defined' ) 1165 #endif 1166 1167 CASE('pno3') 1168 #if defined key_hadocc 1169 ! Dissolved inorganic nitrogen from HadOCC 1170 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1171 #elif defined key_medusa 1172 ! Dissolved inorganic nitrogen from MEDUSA 1173 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1174 #elif defined key_fabm 1175 ! Nitrate from ERSEM 1176 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 1177 #else 1178 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1179 & ' but no biogeochemical model appears to have been defined' ) 1180 #endif 1181 1182 CASE('psi4') 1183 #if defined key_hadocc 1184 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1185 & ' but HadOCC does not simulate silicate' ) 1186 #elif defined key_medusa 1187 ! Silicate from MEDUSA 1188 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1189 #elif defined key_fabm 1190 ! Silicate from ERSEM 1191 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 1192 #else 1193 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1194 & ' but no biogeochemical model appears to have been defined' ) 1195 #endif 1196 1197 CASE('ppo4') 1198 #if defined key_hadocc 1199 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1200 & ' but HadOCC does not simulate phosphate' ) 1201 #elif defined key_medusa 1202 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1203 & ' but MEDUSA does not simulate phosphate' ) 1204 #elif defined key_fabm 1205 ! Phosphate from ERSEM 1206 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 1207 #else 1208 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1209 & ' but no biogeochemical model appears to have been defined' ) 1210 #endif 1211 1212 CASE('pdic') 1213 #if defined key_hadocc 1214 ! Dissolved inorganic carbon from HadOCC 1215 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1216 #elif defined key_medusa 1217 ! Dissolved inorganic carbon from MEDUSA 1218 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1219 #elif defined key_fabm 1220 ! Dissolved inorganic carbon from ERSEM 1221 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 1222 #else 1223 CALL ctl_stop( ' Trying to run pdic observation operator', & 1224 & ' but no biogeochemical model appears to have been defined' ) 1225 #endif 1226 1227 CASE('palk') 1228 #if defined key_hadocc 1229 ! Alkalinity from HadOCC 1230 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1231 #elif defined key_medusa 1232 ! Alkalinity from MEDUSA 1233 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1234 #elif defined key_fabm 1235 ! Alkalinity from ERSEM 1236 zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 1237 #else 1238 CALL ctl_stop( ' Trying to run palk observation operator', & 1239 & ' but no biogeochemical model appears to have been defined' ) 1240 #endif 1241 1242 CASE('pph') 1243 #if defined key_hadocc 1244 CALL ctl_stop( ' Trying to run pph observation operator', & 1245 & ' but HadOCC has no pH diagnostic defined' ) 1246 #elif defined key_medusa && defined key_roam 1247 ! pH from MEDUSA 1248 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1249 #elif defined key_fabm 1250 ! pH from ERSEM 1251 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 1252 #else 1253 CALL ctl_stop( ' Trying to run pph observation operator', & 1254 & ' but no biogeochemical model appears to have been defined' ) 1255 #endif 1256 1257 CASE('po2') 1258 #if defined key_hadocc 1259 CALL ctl_stop( ' Trying to run po2 observation operator', & 1260 & ' but HadOCC does not simulate oxygen' ) 1261 #elif defined key_medusa 1262 ! Oxygen from MEDUSA 1263 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1264 #elif defined key_fabm 1265 ! Oxygen from ERSEM 1266 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 1267 #else 1268 CALL ctl_stop( ' Trying to run po2 observation operator', & 1269 & ' but no biogeochemical model appears to have been defined' ) 1270 #endif 1271 1272 CASE DEFAULT 1273 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1274 1275 END SELECT 1276 1277 DO jvar = 1, profdataqc(jtype)%nvar 1278 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 1279 & nit000, idaystp, jvar, & 1280 & zprofvar(:,:,:,jvar), & 1281 & zprofclim(:,:,:,jvar), & 1282 & fsdept(:,:,:), fsdepw(:,:,:), & 1283 & zprofmask(:,:,:,jvar), & 1284 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1285 & nn_1dint, nn_2dint_default, & 1286 & kdailyavtypes = nn_profdavtypes ) 1287 END DO 1288 1289 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1290 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1291 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1292 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1293 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim ) 1294 1295 END DO 1296 1297 ENDIF 1298 1299 IF ( nsurftypes > 0 ) THEN 1300 1301 !Allocate local work arrays 1302 CALL wrk_alloc( jpi, jpj, zsurfvar ) 1303 CALL wrk_alloc( jpi, jpj, zsurfclim ) 1304 CALL wrk_alloc( jpi, jpj, zsurfmask ) 1305 #if defined key_fabm 1306 CALL wrk_alloc( jpi, jpj, jpk, fabm_3d ) 1307 #endif 1308 1309 DO jtype = 1, nsurftypes 1310 1311 !Defaults which might be changed 1312 zsurfmask(:,:) = tmask(:,:,1) 1313 zsurfclim(:,:) = 0._wp 1314 llog10 = .FALSE. 1315 1316 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 1317 CASE('sst') 1318 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 1319 IF ( ln_output_clim ) zsurfclim(:,:) = tclim(:,:,1) 1320 CASE('sla') 1321 zsurfvar(:,:) = sshn(:,:) 1322 CASE('sss') 1323 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 1324 IF ( ln_output_clim ) zsurfclim(:,:) = sclim(:,:,1) 1325 CASE('sic') 1326 IF ( kstp == 0 ) THEN 1327 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1328 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 1329 & 'time-step but some obs are valid then.' ) 1330 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1331 & ' sea-ice obs will be missed' 1332 ENDIF 1333 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 1334 & surfdataqc(jtype)%nsstp(1) 1335 CYCLE 1336 ELSE 1337 #if defined key_cice 1338 zsurfvar(:,:) = fr_i(:,:) 1339 #elif defined key_lim2 || defined key_lim3 1340 zsurfvar(:,:) = 1._wp - frld(:,:) 1341 #else 1342 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 1343 & ' but no sea-ice model appears to have been defined' ) 1344 #endif 1345 ENDIF 1346 1347 CASE('slchltot') 1348 #if defined key_hadocc 1349 ! Surface chlorophyll from HadOCC 1350 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1351 #elif defined key_medusa 1352 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1353 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1354 #elif defined key_fabm 1355 ! Add all surface chlorophyll groups from ERSEM 1356 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1357 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1358 #else 1359 CALL ctl_stop( ' Trying to run slchltot observation operator', & 1360 & ' but no biogeochemical model appears to have been defined' ) 1361 #endif 1362 llog10 = .TRUE. 1363 1364 CASE('slchldia') 1365 #if defined key_hadocc 1366 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1367 & ' but HadOCC does not explicitly simulate diatoms' ) 1368 #elif defined key_medusa 1369 ! Diatom surface chlorophyll from MEDUSA 1370 zsurfvar(:,:) = trn(:,:,1,jpchd) 1371 #elif defined key_fabm 1372 ! Diatom surface chlorophyll from ERSEM 1373 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 1374 #else 1375 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1376 & ' but no biogeochemical model appears to have been defined' ) 1377 #endif 1378 llog10 = .TRUE. 1379 1380 CASE('slchlnon') 1381 #if defined key_hadocc 1382 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1383 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1384 #elif defined key_medusa 1385 ! Non-diatom surface chlorophyll from MEDUSA 1386 zsurfvar(:,:) = trn(:,:,1,jpchn) 1387 #elif defined key_fabm 1388 ! Add all non-diatom surface chlorophyll groups from ERSEM 1389 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1390 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1391 #else 1392 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1393 & ' but no biogeochemical model appears to have been defined' ) 1394 #endif 1395 llog10 = .TRUE. 1396 1397 CASE('slchldin') 1398 #if defined key_hadocc 1399 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1400 & ' but HadOCC does not explicitly simulate dinoflagellates' ) 1401 #elif defined key_medusa 1402 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1403 & ' but MEDUSA does not explicitly simulate dinoflagellates' ) 1404 #elif defined key_fabm 1405 ! Dinoflagellate surface chlorophyll from ERSEM 1406 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1407 #else 1408 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1409 & ' but no biogeochemical model appears to have been defined' ) 1410 #endif 1411 llog10 = .TRUE. 1412 1413 CASE('slchlmic') 1414 #if defined key_hadocc 1415 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1416 & ' but HadOCC does not explicitly simulate microphytoplankton' ) 1417 #elif defined key_medusa 1418 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1419 & ' but MEDUSA does not explicitly simulate microphytoplankton' ) 1420 #elif defined key_fabm 1421 ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 1422 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1423 #else 1424 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1425 & ' but no biogeochemical model appears to have been defined' ) 1426 #endif 1427 llog10 = .TRUE. 1428 1429 CASE('slchlnan') 1430 #if defined key_hadocc 1431 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1432 & ' but HadOCC does not explicitly simulate nanophytoplankton' ) 1433 #elif defined key_medusa 1434 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1435 & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 1436 #elif defined key_fabm 1437 ! Nanophytoplankton surface chlorophyll from ERSEM 1438 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 1439 #else 1440 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1441 & ' but no biogeochemical model appears to have been defined' ) 1442 #endif 1443 llog10 = .TRUE. 1444 1445 CASE('slchlpic') 1446 #if defined key_hadocc 1447 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1448 & ' but HadOCC does not explicitly simulate picophytoplankton' ) 1449 #elif defined key_medusa 1450 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1451 & ' but MEDUSA does not explicitly simulate picophytoplankton' ) 1452 #elif defined key_fabm 1453 ! Picophytoplankton surface chlorophyll from ERSEM 1454 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 1455 #else 1456 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1457 & ' but no biogeochemical model appears to have been defined' ) 1458 #endif 1459 llog10 = .TRUE. 1460 1461 CASE('schltot') 1462 #if defined key_hadocc 1463 ! Surface chlorophyll from HadOCC 1464 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1465 #elif defined key_medusa 1466 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1467 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1468 #elif defined key_fabm 1469 ! Add all surface chlorophyll groups from ERSEM 1470 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1471 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1472 #else 1473 CALL ctl_stop( ' Trying to run schltot observation operator', & 1474 & ' but no biogeochemical model appears to have been defined' ) 1475 #endif 1476 1477 CASE('slphytot') 1478 #if defined key_hadocc 1479 ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 1480 zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 1481 #elif defined key_medusa 1482 ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 1483 ! multiplied by C:N ratio for each 1484 zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 1485 #elif defined key_fabm 1486 ! Add all surface phytoplankton carbon groups from ERSEM 1487 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1488 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1489 #else 1490 CALL ctl_stop( ' Trying to run slphytot observation operator', & 1491 & ' but no biogeochemical model appears to have been defined' ) 1492 #endif 1493 llog10 = .TRUE. 1494 1495 CASE('slphydia') 1496 #if defined key_hadocc 1497 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1498 & ' but HadOCC does not explicitly simulate diatoms' ) 1499 #elif defined key_medusa 1500 ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1501 zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 1502 #elif defined key_fabm 1503 ! Diatom surface phytoplankton carbon from ERSEM 1504 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 1505 #else 1506 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1507 & ' but no biogeochemical model appears to have been defined' ) 1508 #endif 1509 llog10 = .TRUE. 1510 1511 CASE('slphynon') 1512 #if defined key_hadocc 1513 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1514 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1515 #elif defined key_medusa 1516 ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1517 zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 1518 #elif defined key_fabm 1519 ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 1520 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1521 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1522 #else 1523 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1524 & ' but no biogeochemical model appears to have been defined' ) 1525 #endif 1526 llog10 = .TRUE. 1527 1528 CASE('sspm') 1529 #if defined key_spm 1530 zsurfvar(:,:) = 0.0 1531 DO jn = 1, jp_spm 1532 zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes 1533 END DO 1534 #else 1535 CALL ctl_stop( ' Trying to run sspm observation operator', & 1536 & ' but no spm model appears to have been defined' ) 1537 #endif 1538 1539 CASE('skd490') 1540 #if defined key_hadocc 1541 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1542 & ' but HadOCC does not explicitly simulate Kd490' ) 1543 #elif defined key_medusa 1544 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1545 & ' but MEDUSA does not explicitly simulate Kd490' ) 1546 #elif defined key_fabm 1547 ! light_xEPS diagnostic variable 1548 fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_xeps) 1549 zsurfvar(:,:) = fabm_3d(:,:,1) 1550 #else 1551 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1552 & ' but no biogeochemical model appears to have been defined' ) 1553 #endif 1554 1555 CASE('sfco2') 1556 #if defined key_hadocc 1557 zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 1558 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 1559 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 1560 zsurfvar(:,:) = obfillflt 1561 zsurfmask(:,:) = 0 1562 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 1563 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 1564 ENDIF 1565 #elif defined key_medusa && defined key_roam 1566 zsurfvar(:,:) = f2_fco2w(:,:) 1567 #elif defined key_fabm 1568 ! First, get pCO2 from FABM 1569 fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1570 zsurfvar(:,:) = fabm_3d(:,:,1) 1571 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 1572 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 1573 ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 1574 ! and 1575 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 1576 ! Marine Chemistry, 2: 203-215. 1577 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 1578 ! not explicitly included - atmospheric pressure is not necessarily available so this is 1579 ! the best assumption. 1580 ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 1581 ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 1582 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 1583 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 1584 zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & 1585 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & 1586 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1587 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1588 & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 1589 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 1590 #else 1591 CALL ctl_stop( ' Trying to run sfco2 observation operator', & 1592 & ' but no biogeochemical model appears to have been defined' ) 1593 #endif 1594 1595 CASE('spco2') 1596 #if defined key_hadocc 1597 zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 1598 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 1599 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 1600 zsurfvar(:,:) = obfillflt 1601 zsurfmask(:,:) = 0 1602 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 1603 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 1604 ENDIF 1605 #elif defined key_medusa && defined key_roam 1606 zsurfvar(:,:) = f2_pco2w(:,:) 1607 #elif defined key_fabm 1608 fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1609 zsurfvar(:,:) = fabm_3d(:,:,1) 1610 #else 1611 CALL ctl_stop( ' Trying to run spco2 observation operator', & 1612 & ' but no biogeochemical model appears to have been defined' ) 1613 #endif 1614 1615 CASE DEFAULT 1616 1617 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 1618 1619 END SELECT 1620 1621 IF ( llog10 ) THEN 1622 ! Take the log10 where we can, otherwise exclude 1623 tiny = 1.0e-20 1624 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 1625 zsurfvar(:,:) = LOG10(zsurfvar(:,:)) 1626 ELSEWHERE 1627 zsurfvar(:,:) = obfillflt 1628 zsurfmask(:,:) = 0 1629 END WHERE 1630 ENDIF 1631 1632 IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND. & 1633 & ln_time_mean_sla_bkg ) THEN 1634 !Number of time-steps in meaning period 1635 imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 1636 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1637 & nit000, idaystp, zsurfvar, & 1638 & zsurfclim, zsurfmask, & 1639 & n2dintsurf(jtype), llnightav(jtype), & 1640 & ravglamscl(jtype), ravgphiscl(jtype), & 1641 & lfpindegs(jtype), kmeanstp = imeanstp ) 1642 1061 1643 ELSE 1062 CALL obs_pro_opt( prodatqc(jprofset), & 1063 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1064 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1065 & gdept_1d, tmask, n1dint, n2dint ) 1644 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1645 & nit000, idaystp, zsurfvar, & 1646 & zsurfclim, zsurfmask, & 1647 & n2dintsurf(jtype), llnightav(jtype), & 1648 & ravglamscl(jtype), ravgphiscl(jtype), & 1649 & lfpindegs(jtype) ) 1066 1650 ENDIF 1651 1067 1652 END DO 1653 1654 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 1655 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 1656 #if defined key_fabm 1657 CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d ) 1658 #endif 1659 1068 1660 ENDIF 1069 1661 1070 ! - Sea surface anomaly1071 IF ( ln_sla ) THEN1072 DO jslaset = 1, nslasets1073 CALL obs_sla_opt( sladatqc(jslaset), &1074 & kstp, jpi, jpj, nit000, sshn, &1075 & tmask(:,:,1), n2dint )1076 END DO1077 ENDIF1078 1079 ! - Sea surface temperature1080 IF ( ln_sst ) THEN1081 DO jsstset = 1, nsstsets1082 CALL obs_sst_opt( sstdatqc(jsstset), &1083 & kstp, jpi, jpj, nit000, idaystp, &1084 & tsn(:,:,1,jp_tem), tmask(:,:,1), &1085 & n2dint, ld_sstnight(jsstset) )1086 END DO1087 ENDIF1088 1089 ! - Sea surface salinity1090 IF ( ln_sss ) THEN1091 IF(lwp) WRITE(numout,*) ' SSS currently not available'1092 ENDIF1093 1094 #if defined key_lim2 || defined key_lim31095 IF ( ln_seaice ) THEN1096 DO jseaiceset = 1, nseaicesets1097 CALL obs_seaice_opt( seaicedatqc(jseaiceset), &1098 & kstp, jpi, jpj, nit000, 1.-frld, &1099 & tmask(:,:,1), n2dint )1100 END DO1101 ENDIF1102 #endif1103 1104 ! - Velocity profiles1105 IF ( ln_vel3d ) THEN1106 DO jveloset = 1, nvelosets1107 ! zonal component of velocity1108 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &1109 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, &1110 n1dint, n2dint, ld_velav(jveloset) )1111 END DO1112 ENDIF1113 1114 #if ! defined key_lim2 && ! defined key_lim31115 CALL wrk_dealloc(jpi,jpj,frld)1116 #endif1117 1118 1662 END SUBROUTINE dia_obs 1119 1120 SUBROUTINE dia_obs_wri 1663 1664 SUBROUTINE dia_obs_wri 1121 1665 !!---------------------------------------------------------------------- 1122 1666 !! *** ROUTINE dia_obs_wri *** … … 1126 1670 !! ** Method : Call observation diagnostic output routines 1127 1671 !! 1128 !! ** Action : 1672 !! ** Action : 1129 1673 !! 1130 1674 !! History : … … 1134 1678 !! ! 07-03 (K. Mogensen) General handling of profiles 1135 1679 !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles 1680 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 1136 1681 !!---------------------------------------------------------------------- 1682 !! * Modules used 1683 USE obs_rot_vel ! Rotation of velocities 1684 1137 1685 IMPLICIT NONE 1138 1686 1139 1687 !! * Local declarations 1140 1141 INTEGER :: jprofset ! Profile data set loop variable 1142 INTEGER :: jveloset ! Velocity data set loop variable 1143 INTEGER :: jslaset ! SLA data set loop variable 1144 INTEGER :: jsstset ! SST data set loop variable 1145 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1146 INTEGER :: jset 1147 INTEGER :: jfbini 1148 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1149 CHARACTER(LEN=10) :: cdtmp 1688 INTEGER :: jtype ! Data set loop variable 1689 INTEGER :: jo, jvar, jk 1690 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 1691 & zu, & 1692 & zv 1693 1150 1694 !----------------------------------------------------------------------- 1151 1695 ! Depending on switches call various observation output routines 1152 1696 !----------------------------------------------------------------------- 1153 1697 1154 ! - Temperature/salinity profiles 1155 1156 IF( ln_t3d .OR. ln_s3d ) THEN 1157 1158 ! Copy data from prodatqc to profdata structures 1159 DO jprofset = 1, nprofsets 1160 1161 CALL obs_prof_decompress( prodatqc(jprofset), & 1162 & profdata(jprofset), .TRUE., numout ) 1698 IF ( nproftypes > 0 ) THEN 1699 1700 DO jtype = 1, nproftypes 1701 1702 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 1703 1704 ! For velocity data, rotate the model velocities to N/S, E/W 1705 ! using the compressed data structure. 1706 ALLOCATE( & 1707 & zu(profdataqc(jtype)%nvprot(1)), & 1708 & zv(profdataqc(jtype)%nvprot(2)) & 1709 & ) 1710 1711 CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 1712 1713 DO jo = 1, profdataqc(jtype)%nprof 1714 DO jvar = 1, 2 1715 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 1716 1717 IF ( jvar == 1 ) THEN 1718 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 1719 ELSE 1720 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 1721 ENDIF 1722 1723 END DO 1724 END DO 1725 END DO 1726 1727 DEALLOCATE( zu ) 1728 DEALLOCATE( zv ) 1729 1730 END IF 1731 1732 CALL obs_prof_decompress( profdataqc(jtype), & 1733 & profdata(jtype), .TRUE., numout ) 1734 1735 CALL obs_wri_prof( profdata(jtype) ) 1163 1736 1164 1737 END DO 1165 1738 1166 ! Write the profiles.1167 1168 jprofset = 01169 1170 ! ENACT insitu data1171 1172 IF ( ln_ena ) THEN1173 1174 jprofset = jprofset + 11175 1176 CALL obs_wri_p3d( 'enact', profdata(jprofset) )1177 1178 ENDIF1179 1180 ! Coriolis insitu data1181 1182 IF ( ln_cor ) THEN1183 1184 jprofset = jprofset + 11185 1186 CALL obs_wri_p3d( 'corio', profdata(jprofset) )1187 1188 ENDIF1189 1190 ! Feedback insitu data1191 1192 IF ( ln_profb ) THEN1193 1194 jfbini = jprofset + 11195 1196 DO jprofset = jfbini, nprofsets1197 1198 jset = jprofset - jfbini + 11199 WRITE(cdtmp,'(A,I2.2)')'profb_',jset1200 CALL obs_wri_p3d( cdtmp, profdata(jprofset) )1201 1202 END DO1203 1204 ENDIF1205 1206 1739 ENDIF 1207 1740 1208 ! - Sea surface anomaly1209 IF ( ln_sla ) THEN 1210 1211 ! Copy data from sladatqc to sladata structures 1212 DO jslaset = 1, nslasets1213 1214 CALL obs_surf_decompress( sladatqc(jslaset), & 1215 & sladata(jslaset), .TRUE., numout)1741 IF ( nsurftypes > 0 ) THEN 1742 1743 DO jtype = 1, nsurftypes 1744 1745 CALL obs_surf_decompress( surfdataqc(jtype), & 1746 & surfdata(jtype), .TRUE., numout ) 1747 1748 CALL obs_wri_surf( surfdata(jtype) ) 1216 1749 1217 1750 END DO 1218 1751 1219 jslaset = 01220 1221 ! Write the AVISO SLA data1222 1223 IF ( ln_sladt ) THEN1224 1225 jslaset = 11226 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )1227 jslaset = 21228 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )1229 1230 ENDIF1231 1232 IF ( ln_slafb ) THEN1233 1234 jfbini = jslaset + 11235 1236 DO jslaset = jfbini, nslasets1237 1238 jset = jslaset - jfbini + 11239 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset1240 CALL obs_wri_sla( cdtmp, sladata(jslaset) )1241 1242 END DO1243 1244 ENDIF1245 1246 ENDIF1247 1248 ! - Sea surface temperature1249 IF ( ln_sst ) THEN1250 1251 ! Copy data from sstdatqc to sstdata structures1252 DO jsstset = 1, nsstsets1253 1254 CALL obs_surf_decompress( sstdatqc(jsstset), &1255 & sstdata(jsstset), .TRUE., numout )1256 1257 END DO1258 1259 jsstset = 01260 1261 ! Write the AVISO SST data1262 1263 IF ( ln_reysst ) THEN1264 1265 jsstset = jsstset + 11266 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )1267 1268 ENDIF1269 1270 IF ( ln_ghrsst ) THEN1271 1272 jsstset = jsstset + 11273 CALL obs_wri_sst( 'ghr', sstdata(jsstset) )1274 1275 ENDIF1276 1277 IF ( ln_sstfb ) THEN1278 1279 jfbini = jsstset + 11280 1281 DO jsstset = jfbini, nsstsets1282 1283 jset = jsstset - jfbini + 11284 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset1285 CALL obs_wri_sst( cdtmp, sstdata(jsstset) )1286 1287 END DO1288 1289 ENDIF1290 1291 ENDIF1292 1293 ! - Sea surface salinity1294 IF ( ln_sss ) THEN1295 IF(lwp) WRITE(numout,*) ' SSS currently not available'1296 ENDIF1297 1298 ! - Sea Ice Concentration1299 IF ( ln_seaice ) THEN1300 1301 ! Copy data from seaicedatqc to seaicedata structures1302 DO jseaiceset = 1, nseaicesets1303 1304 CALL obs_surf_decompress( seaicedatqc(jseaiceset), &1305 & seaicedata(jseaiceset), .TRUE., numout )1306 1307 END DO1308 1309 ! Write the Sea Ice data1310 DO jseaiceset = 1, nseaicesets1311 1312 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset1313 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )1314 1315 END DO1316 1317 ENDIF1318 1319 ! Velocity data1320 IF( ln_vel3d ) THEN1321 1322 ! Copy data from veldatqc to velodata structures1323 DO jveloset = 1, nvelosets1324 1325 CALL obs_prof_decompress( veldatqc(jveloset), &1326 & velodata(jveloset), .TRUE., numout )1327 1328 END DO1329 1330 ! Write the profiles.1331 1332 jveloset = 01333 1334 ! Daily averaged data1335 1336 IF ( ln_velavcur ) THEN1337 1338 jveloset = jveloset + 11339 1340 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )1341 1342 ENDIF1343 1344 ! High frequency data1345 1346 IF ( ln_velhrcur ) THEN1347 1348 jveloset = jveloset + 11349 1350 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )1351 1352 ENDIF1353 1354 ! Daily averaged data1355 1356 IF ( ln_velavadcp ) THEN1357 1358 jveloset = jveloset + 11359 1360 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )1361 1362 ENDIF1363 1364 ! High frequency data1365 1366 IF ( ln_velhradcp ) THEN1367 1368 jveloset = jveloset + 11369 1370 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )1371 1372 ENDIF1373 1374 ! Feedback velocity data1375 1376 IF ( ln_velfb ) THEN1377 1378 jfbini = jveloset + 11379 1380 DO jveloset = jfbini, nvelosets1381 1382 jset = jveloset - jfbini + 11383 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset1384 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )1385 1386 END DO1387 1388 ENDIF1389 1390 1752 ENDIF 1391 1753 … … 1405 1767 !! 1406 1768 !!---------------------------------------------------------------------- 1407 ! !obs_grid deallocation1769 ! obs_grid deallocation 1408 1770 CALL obs_grid_deallocate 1409 1771 1410 !! diaobs deallocation 1411 IF ( nprofsets > 0 ) THEN 1412 DEALLOCATE(ld_enact, & 1413 & profdata, & 1414 & prodatqc) 1415 END IF 1416 IF ( ln_sla ) THEN 1417 DEALLOCATE(sladata, & 1418 & sladatqc) 1419 END IF 1420 IF ( ln_seaice ) THEN 1421 DEALLOCATE(sladata, & 1422 & sladatqc) 1423 END IF 1424 IF ( ln_sst ) THEN 1425 DEALLOCATE(sstdata, & 1426 & sstdatqc) 1427 END IF 1428 IF ( ln_vel3d ) THEN 1429 DEALLOCATE(ld_velav, & 1430 & velodata, & 1431 & veldatqc) 1432 END IF 1772 ! diaobs deallocation 1773 IF ( nproftypes > 0 ) & 1774 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 1775 1776 IF ( nsurftypes > 0 ) & 1777 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 1778 & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 1779 1433 1780 END SUBROUTINE dia_obs_dealloc 1434 1781 … … 1436 1783 !!---------------------------------------------------------------------- 1437 1784 !! *** ROUTINE ini_date *** 1438 !!1439 !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format1440 1785 !! 1441 !! ** Method : Get initial datain double precision YYYYMMDD.HHMMSS format1786 !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 1442 1787 !! 1443 !! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format 1788 !! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format 1789 !! 1790 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1444 1791 !! 1445 1792 !! History : … … 1452 1799 USE phycst, ONLY : & ! Physical constants 1453 1800 & rday 1454 ! USE daymod, ONLY : & ! Time variables1455 ! & nmonth_len1456 1801 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1457 1802 & rdt … … 1460 1805 1461 1806 !! * Arguments 1462 REAL( KIND=dp), INTENT(OUT) :: ddobsini! Initial date in YYYYMMDD.HHMMSS1807 REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 1463 1808 1464 1809 !! * Local declarations … … 1468 1813 INTEGER :: ihou 1469 1814 INTEGER :: imin 1470 INTEGER :: imday 1471 REAL(KIND=wp) :: zdayfrc ! Fraction of day1472 1473 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year1474 1475 ! !----------------------------------------------------------------------1476 ! !Initial date initialization (year, month, day, hour, minute)1477 ! !(This assumes that the initial date is for 00z))1478 ! !----------------------------------------------------------------------1815 INTEGER :: imday ! Number of days in month. 1816 INTEGER, DIMENSION(12) :: & 1817 & imonth_len ! Length in days of the months of the current year 1818 REAL(wp) :: zdayfrc ! Fraction of day 1819 1820 !---------------------------------------------------------------------- 1821 ! Initial date initialization (year, month, day, hour, minute) 1822 ! (This assumes that the initial date is for 00z)) 1823 !---------------------------------------------------------------------- 1479 1824 iyea = ndate0 / 10000 1480 1825 imon = ( ndate0 - iyea * 10000 ) / 100 … … 1483 1828 imin = 0 1484 1829 1485 ! !----------------------------------------------------------------------1486 ! !Compute number of days + number of hours + min since initial time1487 ! !----------------------------------------------------------------------1830 !---------------------------------------------------------------------- 1831 ! Compute number of days + number of hours + min since initial time 1832 !---------------------------------------------------------------------- 1488 1833 iday = iday + ( nit000 -1 ) * rdt / rday 1489 1834 zdayfrc = ( nit000 -1 ) * rdt / rday … … 1492 1837 imin = int( (zdayfrc * 24 - ihou) * 60 ) 1493 1838 1494 ! !-----------------------------------------------------------------------1495 ! !Convert number of days (iday) into a real date1496 ! !----------------------------------------------------------------------1839 !----------------------------------------------------------------------- 1840 ! Convert number of days (iday) into a real date 1841 !---------------------------------------------------------------------- 1497 1842 1498 1843 CALL calc_month_len( iyea, imonth_len ) 1499 1844 1500 1845 DO WHILE ( iday > imonth_len(imon) ) 1501 1846 iday = iday - imonth_len(imon) … … 1508 1853 END DO 1509 1854 1510 ! !----------------------------------------------------------------------1511 ! !Convert it into YYYYMMDD.HHMMSS format.1512 ! !----------------------------------------------------------------------1855 !---------------------------------------------------------------------- 1856 ! Convert it into YYYYMMDD.HHMMSS format. 1857 !---------------------------------------------------------------------- 1513 1858 ddobsini = iyea * 10000_dp + imon * 100_dp + & 1514 1859 & iday + ihou * 0.01_dp + imin * 0.0001_dp … … 1520 1865 !!---------------------------------------------------------------------- 1521 1866 !! *** ROUTINE fin_date *** 1522 !!1523 !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format1524 1867 !! 1525 !! ** Method : Get final datain double precision YYYYMMDD.HHMMSS format1868 !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 1526 1869 !! 1527 !! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format 1870 !! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format 1871 !! 1872 !! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format 1528 1873 !! 1529 1874 !! History : … … 1535 1880 USE phycst, ONLY : & ! Physical constants 1536 1881 & rday 1537 ! USE daymod, ONLY : & ! Time variables1538 ! & nmonth_len1539 1882 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1540 1883 & rdt … … 1543 1886 1544 1887 !! * Arguments 1545 REAL( KIND=dp), INTENT(OUT) :: ddobsfin! Final date in YYYYMMDD.HHMMSS1888 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1546 1889 1547 1890 !! * Local declarations … … 1551 1894 INTEGER :: ihou 1552 1895 INTEGER :: imin 1553 INTEGER :: imday 1554 REAL(KIND=wp) :: zdayfrc ! Fraction of day1555 1556 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year1557 1896 INTEGER :: imday ! Number of days in month. 1897 INTEGER, DIMENSION(12) :: & 1898 & imonth_len ! Length in days of the months of the current year 1899 REAL(wp) :: zdayfrc ! Fraction of day 1900 1558 1901 !----------------------------------------------------------------------- 1559 1902 ! Initial date initialization (year, month, day, hour, minute) … … 1565 1908 ihou = 0 1566 1909 imin = 0 1567 1910 1568 1911 !----------------------------------------------------------------------- 1569 1912 ! Compute number of days + number of hours + min since initial time … … 1580 1923 1581 1924 CALL calc_month_len( iyea, imonth_len ) 1582 1925 1583 1926 DO WHILE ( iday > imonth_len(imon) ) 1584 1927 iday = iday - imonth_len(imon) … … 1598 1941 1599 1942 END SUBROUTINE fin_date 1600 1943 1944 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 1945 1946 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1947 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 1948 INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 1949 & ifiles ! Out number of files for each type 1950 CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 1951 & cobstypes ! List of obs types 1952 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 1953 & cfiles ! List of files for all types 1954 1955 !Local variables 1956 INTEGER :: jfile 1957 INTEGER :: jtype 1958 1959 DO jtype = 1, ntypes 1960 1961 ifiles(jtype) = 0 1962 DO jfile = 1, jpmaxnfiles 1963 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 1964 ifiles(jtype) = ifiles(jtype) + 1 1965 END DO 1966 1967 IF ( ifiles(jtype) == 0 ) THEN 1968 CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// & 1969 & ' set to true but no files available to read' ) 1970 ENDIF 1971 1972 IF(lwp) THEN 1973 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1974 DO jfile = 1, ifiles(jtype) 1975 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1976 END DO 1977 ENDIF 1978 1979 END DO 1980 1981 END SUBROUTINE obs_settypefiles 1982 1983 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1984 & n2dint_default, n2dint_type, & 1985 & ravglamscl_type, ravgphiscl_type, & 1986 & lfp_indegs_type, lavnight_type, & 1987 & n2dint, ravglamscl, ravgphiscl, & 1988 & lfpindegs, lavnight ) 1989 1990 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1991 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1992 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1993 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1994 REAL(wp), INTENT(IN) :: & 1995 & ravglamscl_type, & !E/W diameter of obs footprint for this type 1996 & ravgphiscl_type !N/S diameter of obs footprint for this type 1997 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1998 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1999 CHARACTER(len=8), INTENT(IN) :: ctypein 2000 2001 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 2002 & n2dint 2003 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 2004 & ravglamscl, ravgphiscl 2005 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 2006 & lfpindegs, lavnight 2007 2008 lavnight(jtype) = lavnight_type 2009 2010 IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 2011 n2dint(jtype) = n2dint_type 2012 ELSE IF ( n2dint_type == -1 ) THEN 2013 n2dint(jtype) = n2dint_default 2014 ELSE 2015 CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 2016 & ' is not available') 2017 ENDIF 2018 2019 ! For averaging observation footprints set options for size of footprint 2020 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 2021 IF ( ravglamscl_type > 0._wp ) THEN 2022 ravglamscl(jtype) = ravglamscl_type 2023 ELSE 2024 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 2025 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) 2026 ENDIF 2027 2028 IF ( ravgphiscl_type > 0._wp ) THEN 2029 ravgphiscl(jtype) = ravgphiscl_type 2030 ELSE 2031 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 2032 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) 2033 ENDIF 2034 2035 lfpindegs(jtype) = lfp_indegs_type 2036 2037 ENDIF 2038 2039 ! Write out info 2040 IF(lwp) THEN 2041 IF ( n2dint(jtype) <= 4 ) THEN 2042 WRITE(numout,*) ' '//TRIM(ctypein)// & 2043 & ' model counterparts will be interpolated horizontally' 2044 ELSE IF ( n2dint(jtype) <= 6 ) THEN 2045 WRITE(numout,*) ' '//TRIM(ctypein)// & 2046 & ' model counterparts will be averaged horizontally' 2047 WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) 2048 WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) 2049 IF ( lfpindegs(jtype) ) THEN 2050 WRITE(numout,*) ' '//' (in degrees)' 2051 ELSE 2052 WRITE(numout,*) ' '//' (in metres)' 2053 ENDIF 2054 ENDIF 2055 ENDIF 2056 2057 END SUBROUTINE obs_setinterpopts 2058 1601 2059 END MODULE diaobs -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90
r2358 r15670 325 325 CALL obs_mpp_max_integer( kobsj, kobs ) 326 326 ELSE 327 CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj,kobs )327 CALL obs_mpp_find_obs_proc( kproc,kobs ) 328 328 ENDIF 329 329 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90
r8058 r15670 52 52 53 53 !! Default values 54 REAL, PUBLIC :: grid_search_res = 0.5! Resolution of grid54 REAL, PUBLIC :: rn_gridsearchres = 0.5 ! Resolution of grid 55 55 INTEGER, PRIVATE :: gsearch_nlons_def ! Num of longitudes 56 56 INTEGER, PRIVATE :: gsearch_nlats_def ! Num of latitudes … … 83 83 LOGICAL, PUBLIC :: ln_grid_global ! Use global distribution of observations 84 84 CHARACTER(LEN=44), PUBLIC :: & 85 & grid_search_file ! file name head for grid search lookup85 & cn_gridsearchfile ! file name head for grid search lookup 86 86 87 87 !!---------------------------------------------------------------------- … … 613 613 CALL obs_mpp_max_integer( kobsj, kobs ) 614 614 ELSE 615 CALL obs_mpp_find_obs_proc( kproc, kobs i, kobsj, kobs)615 CALL obs_mpp_find_obs_proc( kproc, kobs ) 616 616 ENDIF 617 617 … … 690 690 691 691 IF(lwp) WRITE(numout,*) 692 IF(lwp) WRITE(numout,*)'Grid search resolution : ', grid_search_res693 694 gsearch_nlons_def = NINT( 360.0_wp / grid_search_res )695 gsearch_nlats_def = NINT( 180.0_wp / grid_search_res )696 gsearch_lonmin_def = -180.0_wp + 0.5_wp * grid_search_res697 gsearch_latmin_def = -90.0_wp + 0.5_wp * grid_search_res698 gsearch_dlon_def = grid_search_res699 gsearch_dlat_def = grid_search_res692 IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres 693 694 gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres ) 695 gsearch_nlats_def = NINT( 180.0_wp / rn_gridsearchres ) 696 gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres 697 gsearch_latmin_def = -90.0_wp + 0.5_wp * rn_gridsearchres 698 gsearch_dlon_def = rn_gridsearchres 699 gsearch_dlat_def = rn_gridsearchres 700 700 701 701 IF (lwp) THEN … … 710 710 IF ( ln_grid_global ) THEN 711 711 WRITE(cfname, FMT="(A,'_',A)") & 712 & TRIM( grid_search_file), 'global.nc'712 & TRIM(cn_gridsearchfile), 'global.nc' 713 713 ELSE 714 714 WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") & 715 & TRIM( grid_search_file), nproc, jpni, jpnj715 & TRIM(cn_gridsearchfile), nproc, jpni, jpnj 716 716 ENDIF 717 717 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90
r8058 r15670 35 35 CONTAINS 36 36 37 SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &37 SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 38 38 & pval, pgval, kproc ) 39 39 !!---------------------------------------------------------------------- … … 57 57 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 58 58 INTEGER, INTENT(IN) :: kobs ! Local number of observations 59 INTEGER, INTENT(IN) :: kpi ! Number of points in i direction 60 INTEGER, INTENT(IN) :: kpj ! Number of points in j direction 59 61 INTEGER, INTENT(IN) :: kpk ! Number of levels 60 62 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & … … 63 65 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 64 66 & kproc ! Precomputed processor for each i,j,iobs points 65 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&67 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 66 68 & pval ! Local 3D array to extract data from 67 69 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& … … 73 75 IF (PRESENT(kproc)) THEN 74 76 75 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, &77 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 76 78 & kgrdj, pval, pgval, kproc=kproc ) 77 79 78 80 ELSE 79 81 80 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, &82 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 81 83 & kgrdj, pval, pgval ) 82 84 … … 85 87 ELSE 86 88 87 CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &89 CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 88 90 & pval, pgval ) 89 91 … … 92 94 END SUBROUTINE obs_int_comm_3d 93 95 94 SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, k grdi, kgrdj, pval, pgval, &96 SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kpi, kpj, kgrdi, kgrdj, pval, pgval, & 95 97 & kproc ) 96 98 !!---------------------------------------------------------------------- … … 111 113 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 112 114 INTEGER, INTENT(IN) :: kobs ! Local number of observations 115 INTEGER, INTENT(IN) :: kpi ! Number of model grid points in i direction 116 INTEGER, INTENT(IN) :: kpj ! Number of model grid points in j direction 113 117 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 114 118 & kgrdi, & ! i,j indicies for each stencil … … 116 120 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 117 121 & kproc ! Precomputed processor for each i,j,iobs points 118 REAL(KIND=wp), DIMENSION( jpi,jpj), INTENT(IN) ::&122 REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) ::& 119 123 & pval ! Local 3D array to extra data from 120 124 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::& … … 136 140 IF (PRESENT(kproc)) THEN 137 141 138 CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &142 CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 139 143 & zgval, kproc=kproc ) 140 144 ELSE 141 145 142 CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &146 CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 143 147 & zgval ) 144 148 … … 154 158 END SUBROUTINE obs_int_comm_2d 155 159 156 SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &160 SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 157 161 & pval, pgval, kproc ) 158 162 !!---------------------------------------------------------------------- … … 174 178 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 175 179 INTEGER, INTENT(IN) :: kobs ! Local number of observations 180 INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction 181 INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction 176 182 INTEGER, INTENT(IN) :: kpk ! Number of levels 177 183 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & … … 180 186 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 181 187 & kproc ! Precomputed processor for each i,j,iobs points 182 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&188 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 183 189 & pval ! Local 3D array to extract data from 184 190 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& … … 207 213 208 214 ! Check valid points 209 215 210 216 IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. & 211 217 & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN 212 218 213 219 CALL ctl_stop( 'Error in obs_int_comm_3d_global', & 214 220 & 'Point outside global domain' ) 215 221 216 222 ENDIF 217 223 … … 323 329 END SUBROUTINE obs_int_comm_3d_global 324 330 325 SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &331 SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 326 332 & pval, pgval ) 327 333 !!---------------------------------------------------------------------- … … 343 349 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 344 350 INTEGER, INTENT(IN) :: kobs ! Local number of observations 351 INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction 352 INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction 345 353 INTEGER, INTENT(IN) :: kpk ! Number of levels 346 354 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 347 355 & kgrdi, & ! i,j indicies for each stencil 348 356 & kgrdj 349 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&357 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 350 358 & pval ! Local 3D array to extract data from 351 359 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90
r8058 r15670 13 13 !! ** Method : Straightforward search 14 14 !! 15 !! ** Action : 15 !! ** Action : Will return level associated with T-point below the obs 16 !! depth, except when observation is in the top box will 17 !! return level 2. Also, if obs depth greater than depth 18 !! of last wet T-point (kpk-1) will return level kpk. 16 19 !! 17 20 !! History : … … 43 46 DO ji = 1, kobs 44 47 kobsk(ji) = 1 45 depk: DO jk = 2, kgrd 46 IF ( pgrddep(jk) > =pobsdep(ji) ) EXIT depk48 depk: DO jk = 2, kgrd-1 49 IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk 47 50 END DO depk 48 51 kobsk(ji) = jk -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90
r8058 r15670 7 7 !! - ! 2006-05 (K. Mogensen) Reformatted 8 8 !! - ! 2008-01 (K. Mogensen) add mpp_global_max 9 !! 3.6 ! 2015-01 (J. Waters) obs_mpp_find_obs_proc 10 !! rewritten to avoid global arrays 9 11 !!---------------------------------------------------------------------- 10 12 # define mpivar mpi_double_precision … … 12 14 !! obs_mpp_bcast_integer : Broadcast an integer array from a processor to all processors 13 15 !! obs_mpp_max_integer : Find maximum on all processors of each value in an integer on all processors 14 !! obs_mpp_find_obs_proc : Find processors which should hold the observations 16 !! obs_mpp_find_obs_proc : Find processors which should hold the observations, avoiding global arrays 15 17 !! obs_mpp_sum_integers : Sum an integer array from all processors 16 18 !! obs_mpp_sum_integer : Sum an integer from all processors … … 96 98 ! 97 99 INTEGER :: ierr 98 INTEGER, DIMENSION(kno) :: ivals 99 ! 100 INCLUDE 'mpif.h' 101 !!---------------------------------------------------------------------- 100 INTEGER, DIMENSION(:), ALLOCATABLE :: ivals 101 ! 102 INCLUDE 'mpif.h' 103 !!---------------------------------------------------------------------- 104 105 ALLOCATE( ivals(kno) ) 102 106 103 107 ! Call the MPI library to find the maximum across processors … … 105 109 & mpi_max, mpi_comm_opa, ierr ) 106 110 kvals(:) = ivals(:) 111 112 DEALLOCATE( ivals ) 107 113 #else 108 114 ! no MPI: empty routine … … 111 117 112 118 113 SUBROUTINE obs_mpp_find_obs_proc( kobsp, kobsi, kobsj,kno )114 !!---------------------------------------------------------------------- 115 !! *** ROUTINE obs_mpp_find_obs_proc ***116 !! 117 !! ** Purpose : From the array kobsp containing the results of the grid119 SUBROUTINE obs_mpp_find_obs_proc( kobsp,kno ) 120 !!---------------------------------------------------------------------- 121 !! *** ROUTINE obs_mpp_find_obs_proc *** 122 !! 123 !! ** Purpose : From the array kobsp containing the results of the 118 124 !! grid search on each processor the processor return a 119 125 !! decision of which processors should hold the observation. 120 126 !! 121 !! ** Method : A temporary 2D array holding all the decisions is122 !! constructed using mpi_allgather on each processor.123 !! If more than one processor has found the observation124 !! with the observation in the inner domain gets it125 !! 126 !! ** Action : This does only work for MPI. 127 !! ** Method : Synchronize the processor number for each obs using 128 !! obs_mpp_max_integer. If an observation exists on two 129 !! processors it will be allocated to the lower numbered 130 !! processor. 131 !! 132 !! ** Action : This does only work for MPI. 127 133 !! It does not work for SHMEM. 128 134 !! … … 130 136 !!---------------------------------------------------------------------- 131 137 INTEGER , INTENT(in ) :: kno 132 INTEGER, DIMENSION(kno), INTENT(in ) :: kobsi, kobsj133 138 INTEGER, DIMENSION(kno), INTENT(inout) :: kobsp 134 139 ! 135 140 #if defined key_mpp_mpi 136 141 ! 137 INTEGER :: ji 138 INTEGER :: jj 139 INTEGER :: size 140 INTEGER :: ierr 141 INTEGER :: iobsip 142 INTEGER :: iobsjp 143 INTEGER :: num_sus_obs 144 INTEGER, DIMENSION(kno) :: iobsig, iobsjg 145 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iobsp, iobsi, iobsj 146 !! 147 INCLUDE 'mpif.h' 148 !!---------------------------------------------------------------------- 149 150 !----------------------------------------------------------------------- 151 ! Call the MPI library to find the maximum accross processors 152 !----------------------------------------------------------------------- 153 CALL mpi_comm_size( mpi_comm_opa, size, ierr ) 154 !----------------------------------------------------------------------- 155 ! Convert local grids points to global grid points 156 !----------------------------------------------------------------------- 142 ! 143 INTEGER :: ji, isum 144 INTEGER, DIMENSION(:), ALLOCATABLE :: iobsp 145 !! 146 !! 147 148 ALLOCATE( iobsp(kno) ) 149 150 iobsp(:)=kobsp(:) 151 152 WHERE( iobsp(:) == -1 ) 153 iobsp(:) = 9999999 154 END WHERE 155 156 iobsp(:)=-1*iobsp(:) 157 158 CALL obs_mpp_max_integer( iobsp, kno ) 159 160 kobsp(:)=-1*iobsp(:) 161 162 isum=0 157 163 DO ji = 1, kno 158 IF ( ( kobsi(ji) >= 1 ) .AND. ( kobsi(ji) <= jpi ) .AND. & 159 & ( kobsj(ji) >= 1 ) .AND. ( kobsj(ji) <= jpj ) ) THEN 160 iobsig(ji) = mig( kobsi(ji) ) 161 iobsjg(ji) = mjg( kobsj(ji) ) 162 ELSE 163 iobsig(ji) = -1 164 iobsjg(ji) = -1 164 IF ( kobsp(ji) == 9999999 ) THEN 165 isum=isum+1 166 kobsp(ji)=-1 165 167 ENDIF 166 END DO 167 !----------------------------------------------------------------------- 168 ! Get the decisions from all processors 169 !----------------------------------------------------------------------- 170 ALLOCATE( iobsp(kno,size) ) 171 ALLOCATE( iobsi(kno,size) ) 172 ALLOCATE( iobsj(kno,size) ) 173 CALL mpi_allgather( kobsp, kno, mpi_integer, & 174 & iobsp, kno, mpi_integer, & 175 & mpi_comm_opa, ierr ) 176 CALL mpi_allgather( iobsig, kno, mpi_integer, & 177 & iobsi, kno, mpi_integer, & 178 & mpi_comm_opa, ierr ) 179 CALL mpi_allgather( iobsjg, kno, mpi_integer, & 180 & iobsj, kno, mpi_integer, & 181 & mpi_comm_opa, ierr ) 182 183 !----------------------------------------------------------------------- 184 ! Find the processor with observations from the lowest processor 185 ! number among processors holding the observation. 186 !----------------------------------------------------------------------- 187 kobsp(:) = -1 188 num_sus_obs = 0 189 DO ji = 1, kno 190 DO jj = 1, size 191 IF ( ( kobsp(ji) == -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 192 kobsp(ji) = iobsp(ji,jj) 193 iobsip = iobsi(ji,jj) 194 iobsjp = iobsj(ji,jj) 195 ENDIF 196 IF ( ( kobsp(ji) /= -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 197 IF ( ( iobsip /= iobsi(ji,jj) ) .OR. & 198 & ( iobsjp /= iobsj(ji,jj) ) ) THEN 199 IF ( ( kobsp(ji) < 1000000 ) .AND. & 200 & ( iobsp(ji,jj) < 1000000 ) ) THEN 201 num_sus_obs=num_sus_obs+1 202 ENDIF 203 ENDIF 204 IF ( mppmap(iobsip,iobsjp) /= ( kobsp(ji)+1 ) ) THEN 205 IF ( ( iobsi(ji,jj) /= -1 ) .AND. & 206 & ( iobsj(ji,jj) /= -1 ) ) THEN 207 IF ((mppmap(iobsi(ji,jj),iobsj(ji,jj)) == (iobsp(ji,jj)+1))& 208 & .OR. ( iobsp(ji,jj) < kobsp(ji) ) ) THEN 209 kobsp(ji) = iobsp(ji,jj) 210 iobsip = iobsi(ji,jj) 211 iobsjp = iobsj(ji,jj) 212 ENDIF 213 ENDIF 214 ENDIF 215 ENDIF 216 END DO 217 END DO 218 IF (lwp) WRITE(numout,*) 'Number of suspicious observations: ',num_sus_obs 219 220 DEALLOCATE( iobsj ) 221 DEALLOCATE( iobsi ) 168 ENDDO 169 170 171 IF ( isum > 0 ) THEN 172 IF (lwp) WRITE(numout,*) isum, ' observations failed the grid search.' 173 IF (lwp) WRITE(numout,*)'If ln_grid_search_lookup=.TRUE., try reducing grid_search_res' 174 ENDIF 175 222 176 DEALLOCATE( iobsp ) 177 223 178 #else 224 179 ! no MPI: empty routine 225 #endif 226 !180 #endif 181 227 182 END SUBROUTINE obs_mpp_find_obs_proc 228 183 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r8058 r15670 7 7 8 8 !!---------------------------------------------------------------------- 9 !! obs_pro_opt : Compute the model counterpart of temperature and 10 !! salinity observations from profiles 11 !! obs_sla_opt : Compute the model counterpart of sea level anomaly 12 !! observations 13 !! obs_sst_opt : Compute the model counterpart of sea surface temperature 14 !! observations 15 !! obs_sss_opt : Compute the model counterpart of sea surface salinity 16 !! observations 17 !! obs_seaice_opt : Compute the model counterpart of sea ice concentration 18 !! observations 19 !! 20 !! obs_vel_opt : Compute the model counterpart of zonal and meridional 21 !! components of velocity from observations. 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 !! obs_surf_opt : Compute the model counterpart of surface data 22 11 !!---------------------------------------------------------------------- 23 12 24 !! * Modules used 13 !! * Modules used 25 14 USE par_kind, ONLY : & ! Precision variables 26 15 & wp 27 16 USE in_out_manager ! I/O manager 28 17 USE obs_inter_sup ! Interpolation support 29 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs ervationpt18 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 30 19 & obs_int_h2d, & 31 20 & obs_int_h2d_init 32 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 33 26 & obs_int_z1d, & 34 27 & obs_int_z1d_spl 35 USE obs_const, ONLY : &36 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 37 30 USE dom_oce, ONLY : & 38 & glamt, glam u, glamv, &39 & gphit, gphi u, gphiv40 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 41 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 37 USE obs_grid, ONLY : & 38 & obs_level_search 42 39 43 40 IMPLICIT NONE … … 46 43 PRIVATE 47 44 48 PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations 49 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 50 & obs_sst_opt, & ! Compute the model counterpart of SST observations 51 & obs_sss_opt, & ! Compute the model counterpart of SSS observations 52 & obs_seaice_opt, & 53 & obs_vel_opt ! Compute the model counterpart of velocity profile data 54 55 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_surf_opt ! Compute the model counterpart of surface obs 47 48 INTEGER, PARAMETER, PUBLIC :: & 49 & imaxavtypes = 20 ! Max number of daily avgd obs types 56 50 57 51 !!---------------------------------------------------------------------- … … 61 55 !!---------------------------------------------------------------------- 62 56 57 !! * Substitutions 58 # include "domzgr_substitute.h90" 63 59 CONTAINS 64 60 65 SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 66 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 67 & kdailyavtypes ) 61 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 63 & kit000, kdaystp, kvar, & 64 & pvar, pclim, & 65 & pgdept, pgdepw, pmask, & 66 & plam, pphi, & 67 & k1dint, k2dint, kdailyavtypes ) 68 68 69 !!----------------------------------------------------------------------- 69 70 !! … … 78 79 !! 79 80 !! First, a vertical profile of horizontally interpolated model 80 !! now temperatures is computed at the obs (lon, lat) point.81 !! now values is computed at the obs (lon, lat) point. 81 82 !! Several horizontal interpolation schemes are available: 82 83 !! - distance-weighted (great circle) (k2dint = 0) … … 86 87 !! - polynomial (quadrilateral grid) (k2dint = 4) 87 88 !! 88 !! Next, the vertical temperatureprofile is interpolated to the89 !! Next, the vertical profile is interpolated to the 89 90 !! data depth points. Two vertical interpolation schemes are 90 91 !! available: … … 96 97 !! routine. 97 98 !! 98 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is99 !! If the logical is switched on, the model equivalent is 99 100 !! a daily mean model temperature field. So, we first compute 100 101 !! the mean, then interpolate only at the end of the day. 101 102 !! 102 !! Note: thein situ temperature observations must be converted103 !! Note: in situ temperature observations must be converted 103 104 !! to potential temperature (the model variable) prior to 104 105 !! assimilation. 105 !!??????????????????????????????????????????????????????????????106 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???107 !!??????????????????????????????????????????????????????????????108 106 !! 109 107 !! ** Action : … … 115 113 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 116 114 !! ! 07-03 (K. Mogensen) General handling of profiles 115 !! ! 15-02 (M. Martin) Combined routine for all profile types 116 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 117 117 !!----------------------------------------------------------------------- 118 118 119 119 !! * Modules used 120 120 USE obs_profiles_def ! Definition of storage space for profile obs. … … 123 123 124 124 !! * Arguments 125 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 126 INTEGER, INTENT(IN) :: kt ! Time step 127 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 125 TYPE(obs_prof), INTENT(INOUT) :: & 126 & prodatqc ! Subset of profile data passing QC 127 INTEGER, INTENT(IN) :: kt ! Time step 128 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 128 129 INTEGER, INTENT(IN) :: kpj 129 130 INTEGER, INTENT(IN) :: kpk 130 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 131 ! (kit000-1 = restart time) 132 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 133 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 134 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 131 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 132 ! (kit000-1 = restart time) 133 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 134 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 135 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 136 INTEGER, INTENT(IN) :: kvar ! Number of variable in prodatqc 135 137 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 136 & ptn, & ! Model temperature field 137 & psn, & ! Model salinity field 138 & ptmask ! Land-sea mask 139 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 140 & pgdept ! Model array of depth levels 138 & pvar, & ! Model field for variable 139 & pclim, & ! Climatology field for variable 140 & pmask ! Land-sea mask for variable 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 142 & plam, & ! Model longitudes for variable 143 & pphi ! Model latitudes for variable 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 145 & pgdept, & ! Model array of depth T levels 146 & pgdepw ! Model array of depth W levels 141 147 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 142 & kdailyavtypes! Types for daily averages 148 & kdailyavtypes ! Types for daily averages 149 143 150 !! * Local declarations 144 151 INTEGER :: ji … … 152 159 INTEGER :: iend 153 160 INTEGER :: iobs 161 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 162 INTEGER :: inum_obs 154 163 INTEGER, DIMENSION(imaxavtypes) :: & 155 164 & idailyavtypes 165 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 166 & igrdi, & 167 & igrdj 168 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 169 156 170 REAL(KIND=wp) :: zlam 157 171 REAL(KIND=wp) :: zphi 158 172 REAL(KIND=wp) :: zdaystp 159 173 REAL(KIND=wp), DIMENSION(kpk) :: & 160 & zobsmask, &161 174 & zobsk, & 162 & zobs2k 163 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 175 & zobs2k, & 176 & zclm2k 177 REAL(KIND=wp), DIMENSION(2,2,1) :: & 178 & zweig1, & 164 179 & zweig 165 180 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 166 & zmask, & 167 & zintt, & 168 & zints, & 169 & zinmt, & 170 & zinms 181 & zmask, & 182 & zclim, & 183 & zint, & 184 & zinm, & 185 & zgdept, & 186 & zgdepw 171 187 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 172 & zglam, &188 & zglam, & 173 189 & zgphi 174 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 175 & igrdi, & 176 & igrdj 190 REAL(KIND=wp), DIMENSION(1) :: zmsk 191 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 192 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim 193 194 LOGICAL :: ld_dailyav 177 195 178 196 !------------------------------------------------------------------------ 179 197 ! Local initialization 180 198 !------------------------------------------------------------------------ 181 ! ...Record and data counters199 ! Record and data counters 182 200 inrc = kt - kit000 + 2 183 201 ipro = prodatqc%npstp(inrc) 184 202 185 203 ! Daily average types 204 ld_dailyav = .FALSE. 186 205 IF ( PRESENT(kdailyavtypes) ) THEN 187 206 idailyavtypes(:) = kdailyavtypes(:) 207 IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 188 208 ELSE 189 209 idailyavtypes(:) = -1 190 210 ENDIF 191 211 192 ! Initialize daily mean for first timestep 212 ! Daily means are calculated for values over timesteps: 213 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 193 214 idayend = MOD( kt - kit000 + 1, kdaystp ) 194 215 195 ! Added kt == 0 test to catch restart case 196 IF ( idayend == 1 .OR. kt == 0) THEN 197 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 216 IF ( ld_dailyav ) THEN 217 218 ! Initialize daily mean for first timestep of the day 219 IF ( idayend == 1 .OR. kt == 0 ) THEN 220 DO jk = 1, jpk 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 224 END DO 225 END DO 226 END DO 227 ENDIF 228 198 229 DO jk = 1, jpk 199 230 DO jj = 1, jpj 200 231 DO ji = 1, jpi 201 prodatqc%vdmean(ji,jj,jk,1) = 0.0 202 prodatqc%vdmean(ji,jj,jk,2) = 0.0 232 ! Increment field for computing daily mean 233 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 234 & + pvar(ji,jj,jk) 203 235 END DO 204 236 END DO 205 237 END DO 206 ENDIF 207 208 DO jk = 1, jpk 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 ! Increment the temperature field for computing daily mean 212 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 213 & + ptn(ji,jj,jk) 214 ! Increment the salinity field for computing daily mean 215 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 216 & + psn(ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 ! Compute the daily mean at the end of day 222 zdaystp = 1.0 / REAL( kdaystp ) 223 IF ( idayend == 0 ) THEN 224 DO jk = 1, jpk 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 228 & * zdaystp 229 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 230 & * zdaystp 238 239 ! Compute the daily mean at the end of day 240 zdaystp = 1.0 / REAL( kdaystp ) 241 IF ( idayend == 0 ) THEN 242 IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 243 CALL FLUSH(numout) 244 DO jk = 1, jpk 245 DO jj = 1, jpj 246 DO ji = 1, jpi 247 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 248 & * zdaystp 249 END DO 231 250 END DO 232 251 END DO 233 END DO 252 ENDIF 253 234 254 ENDIF 235 255 236 256 ! Get the data for interpolation 237 257 ALLOCATE( & 238 & igrdi(2,2,ipro), & 239 & igrdj(2,2,ipro), & 240 & zglam(2,2,ipro), & 241 & zgphi(2,2,ipro), & 242 & zmask(2,2,kpk,ipro), & 243 & zintt(2,2,kpk,ipro), & 244 & zints(2,2,kpk,ipro) & 258 & igrdi(2,2,ipro), & 259 & igrdj(2,2,ipro), & 260 & zglam(2,2,ipro), & 261 & zgphi(2,2,ipro), & 262 & zmask(2,2,kpk,ipro), & 263 & zint(2,2,kpk,ipro), & 264 & zgdept(2,2,kpk,ipro), & 265 & zgdepw(2,2,kpk,ipro) & 245 266 & ) 267 268 IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) ) 246 269 247 270 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 248 271 iobs = jobs - prodatqc%nprofup 249 igrdi(1,1,iobs) = prodatqc%mi(jobs, 1)-1250 igrdj(1,1,iobs) = prodatqc%mj(jobs, 1)-1251 igrdi(1,2,iobs) = prodatqc%mi(jobs, 1)-1252 igrdj(1,2,iobs) = prodatqc%mj(jobs, 1)253 igrdi(2,1,iobs) = prodatqc%mi(jobs, 1)254 igrdj(2,1,iobs) = prodatqc%mj(jobs, 1)-1255 igrdi(2,2,iobs) = prodatqc%mi(jobs, 1)256 igrdj(2,2,iobs) = prodatqc%mj(jobs, 1)272 igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 273 igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 274 igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 275 igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 276 igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 277 igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 278 igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 279 igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 257 280 END DO 258 281 259 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 260 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 261 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 262 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 263 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 264 282 ! Initialise depth arrays 283 zgdept(:,:,:,:) = 0.0 284 zgdepw(:,:,:,:) = 0.0 285 286 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 287 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 288 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 289 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint ) 290 291 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 292 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 293 294 IF ( prodatqc%lclim ) THEN 295 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim ) 296 ENDIF 297 265 298 ! At the end of the day also get interpolated means 266 IF ( idayend == 0 ) THEN267 268 ALLOCATE( &269 & zinmt(2,2,kpk,ipro), & 270 & zinms(2,2,kpk,ipro)&271 & )272 273 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &274 & prodatqc%vdmean(:,:,:,1), zinmt ) 275 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &276 & prodatqc%vdmean(:,:,:,2), zinms )277 278 ENDIF299 IF ( ld_dailyav .AND. idayend == 0 ) THEN 300 301 ALLOCATE( zinm(2,2,kpk,ipro) ) 302 303 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 304 & prodatqc%vdmean(:,:,:,kvar), zinm ) 305 306 ENDIF 307 308 ! Return if no observations to process 309 ! Has to be done after comm commands to ensure processors 310 ! stay in sync 311 IF ( ipro == 0 ) RETURN 279 312 280 313 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro … … 283 316 284 317 IF ( kt /= prodatqc%mstp(jobs) ) THEN 285 318 286 319 IF(lwp) THEN 287 320 WRITE(numout,*) … … 298 331 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 299 332 ENDIF 300 333 301 334 zlam = prodatqc%rlam(jobs) 302 335 zphi = prodatqc%rphi(jobs) 336 337 ! Horizontal weights 338 ! Masked values are calculated later. 339 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 340 341 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 342 & zglam(:,:,iobs), zgphi(:,:,iobs), & 343 & zmask(:,:,1,iobs), zweig1, zmsk ) 344 345 ENDIF 346 347 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 348 349 zobsk(:) = obfillflt 350 351 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 352 353 IF ( idayend == 0 ) THEN 354 ! Daily averaged data 355 356 ! vertically interpolate all 4 corners 357 ista = prodatqc%npvsta(jobs,kvar) 358 iend = prodatqc%npvend(jobs,kvar) 359 inum_obs = iend - ista + 1 360 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 361 IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 362 363 DO iin=1,2 364 DO ijn=1,2 365 366 IF ( k1dint == 1 ) THEN 367 CALL obs_int_z1d_spl( kpk, & 368 & zinm(iin,ijn,:,iobs), & 369 & zobs2k, zgdept(iin,ijn,:,iobs), & 370 & zmask(iin,ijn,:,iobs)) 371 372 IF ( prodatqc%lclim ) THEN 373 CALL obs_int_z1d_spl( kpk, & 374 & zclim(iin,ijn,:,iobs), & 375 & zclm2k, zgdept(iin,ijn,:,iobs), & 376 & zmask(iin,ijn,:,iobs)) 377 ENDIF 378 379 ENDIF 380 381 CALL obs_level_search(kpk, & 382 & zgdept(iin,ijn,:,iobs), & 383 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 384 & iv_indic) 385 386 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 387 & prodatqc%var(kvar)%vdep(ista:iend), & 388 & zinm(iin,ijn,:,iobs), & 389 & zobs2k, interp_corner(iin,ijn,:), & 390 & zgdept(iin,ijn,:,iobs), & 391 & zmask(iin,ijn,:,iobs)) 392 393 IF ( prodatqc%lclim ) THEN 394 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 395 & prodatqc%var(kvar)%vdep(ista:iend), & 396 & zclim(iin,ijn,:,iobs), & 397 & zclm2k, interp_corner_clim(iin,ijn,:), & 398 & zgdept(iin,ijn,:,iobs), & 399 & zmask(iin,ijn,:,iobs)) 400 ENDIF 401 402 ENDDO 403 ENDDO 404 405 ENDIF !idayend 406 407 ELSE 408 409 ! Point data 410 411 ! vertically interpolate all 4 corners 412 ista = prodatqc%npvsta(jobs,kvar) 413 iend = prodatqc%npvend(jobs,kvar) 414 inum_obs = iend - ista + 1 415 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 416 IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 417 DO iin=1,2 418 DO ijn=1,2 419 420 IF ( k1dint == 1 ) THEN 421 CALL obs_int_z1d_spl( kpk, & 422 & zint(iin,ijn,:,iobs),& 423 & zobs2k, zgdept(iin,ijn,:,iobs), & 424 & zmask(iin,ijn,:,iobs)) 425 426 IF ( prodatqc%lclim ) THEN 427 CALL obs_int_z1d_spl( kpk, & 428 & zclim(iin,ijn,:,iobs),& 429 & zclm2k, zgdept(iin,ijn,:,iobs), & 430 & zmask(iin,ijn,:,iobs)) 431 ENDIF 432 433 ENDIF 434 435 CALL obs_level_search(kpk, & 436 & zgdept(iin,ijn,:,iobs),& 437 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 438 & iv_indic) 439 440 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 441 & prodatqc%var(kvar)%vdep(ista:iend), & 442 & zint(iin,ijn,:,iobs), & 443 & zobs2k,interp_corner(iin,ijn,:), & 444 & zgdept(iin,ijn,:,iobs), & 445 & zmask(iin,ijn,:,iobs) ) 446 447 IF ( prodatqc%lclim ) THEN 448 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 449 & prodatqc%var(kvar)%vdep(ista:iend), & 450 & zclim(iin,ijn,:,iobs), & 451 & zclm2k,interp_corner_clim(iin,ijn,:), & 452 & zgdept(iin,ijn,:,iobs), & 453 & zmask(iin,ijn,:,iobs) ) 454 ENDIF 303 455 304 ! Horizontal weights and vertical mask 305 306 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 307 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 308 309 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 310 & zglam(:,:,iobs), zgphi(:,:,iobs), & 311 & zmask(:,:,:,iobs), zweig, zobsmask ) 312 313 ENDIF 314 315 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 316 317 zobsk(:) = obfillflt 318 319 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 320 321 IF ( idayend == 0 ) THEN 456 ENDDO 457 ENDDO 458 459 ENDIF 460 461 !------------------------------------------------------------- 462 ! Compute the horizontal interpolation for every profile level 463 !------------------------------------------------------------- 464 465 DO ikn=1,inum_obs 466 iend=ista+ikn-1 322 467 323 ! Daily averaged moored buoy (MRB) data 324 325 CALL obs_int_h2d( kpk, kpk, & 326 & zweig, zinmt(:,:,:,iobs), zobsk ) 327 328 329 ELSE 330 331 CALL ctl_stop( ' A nonzero' // & 332 & ' number of profile T BUOY data should' // & 333 & ' only occur at the end of a given day' ) 334 468 zweig(:,:,1) = 0._wp 469 470 ! This code forces the horizontal weights to be 471 ! zero IF the observation is below the bottom of the 472 ! corners of the interpolation nodes, Or if it is in 473 ! the mask. This is important for observations near 474 ! steep bathymetry 475 DO iin=1,2 476 DO ijn=1,2 477 478 depth_loop: DO ik=kpk,2,-1 479 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 480 481 zweig(iin,ijn,1) = & 482 & zweig1(iin,ijn,1) * & 483 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 484 & - prodatqc%var(kvar)%vdep(iend)),0._wp) 485 486 EXIT depth_loop 487 488 ENDIF 489 490 ENDDO depth_loop 491 492 ENDDO 493 ENDDO 494 495 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 496 & prodatqc%var(kvar)%vmod(iend:iend) ) 497 498 IF ( prodatqc%lclim ) THEN 499 CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 500 & prodatqc%var(kvar)%vclm(iend:iend) ) 335 501 ENDIF 336 337 ELSE 338 339 ! Point data 340 341 CALL obs_int_h2d( kpk, kpk, & 342 & zweig, zintt(:,:,:,iobs), zobsk ) 343 344 ENDIF 345 346 !------------------------------------------------------------- 347 ! Compute vertical second-derivative of the interpolating 348 ! polynomial at obs points 349 !------------------------------------------------------------- 350 351 IF ( k1dint == 1 ) THEN 352 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 353 & pgdept, zobsmask ) 354 ENDIF 355 356 !----------------------------------------------------------------- 357 ! Vertical interpolation to the observation point 358 !----------------------------------------------------------------- 359 ista = prodatqc%npvsta(jobs,1) 360 iend = prodatqc%npvend(jobs,1) 361 CALL obs_int_z1d( kpk, & 362 & prodatqc%var(1)%mvk(ista:iend), & 363 & k1dint, iend - ista + 1, & 364 & prodatqc%var(1)%vdep(ista:iend), & 365 & zobsk, zobs2k, & 366 & prodatqc%var(1)%vmod(ista:iend), & 367 & pgdept, zobsmask ) 368 369 ENDIF 370 371 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 372 373 zobsk(:) = obfillflt 374 375 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 376 377 IF ( idayend == 0 ) THEN 378 379 ! Daily averaged moored buoy (MRB) data 380 381 CALL obs_int_h2d( kpk, kpk, & 382 & zweig, zinms(:,:,:,iobs), zobsk ) 383 384 ELSE 385 386 CALL ctl_stop( ' A nonzero' // & 387 & ' number of profile S BUOY data should' // & 388 & ' only occur at the end of a given day' ) 389 390 ENDIF 391 392 ELSE 393 394 ! Point data 395 396 CALL obs_int_h2d( kpk, kpk, & 397 & zweig, zints(:,:,:,iobs), zobsk ) 398 399 ENDIF 400 401 402 !------------------------------------------------------------- 403 ! Compute vertical second-derivative of the interpolating 404 ! polynomial at obs points 405 !------------------------------------------------------------- 406 407 IF ( k1dint == 1 ) THEN 408 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 409 & pgdept, zobsmask ) 410 ENDIF 411 412 !---------------------------------------------------------------- 413 ! Vertical interpolation to the observation point 414 !---------------------------------------------------------------- 415 ista = prodatqc%npvsta(jobs,2) 416 iend = prodatqc%npvend(jobs,2) 417 CALL obs_int_z1d( kpk, & 418 & prodatqc%var(2)%mvk(ista:iend),& 419 & k1dint, iend - ista + 1, & 420 & prodatqc%var(2)%vdep(ista:iend),& 421 & zobsk, zobs2k, & 422 & prodatqc%var(2)%vmod(ista:iend),& 423 & pgdept, zobsmask ) 424 425 ENDIF 426 427 END DO 502 503 ! Set QC flag for any observations found below the bottom 504 ! needed as the check here is more strict than that in obs_prep 505 IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 428 506 507 ENDDO 508 509 DEALLOCATE(interp_corner,iv_indic) 510 IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim ) 511 512 ENDIF 513 514 ENDDO 515 429 516 ! Deallocate the data for interpolation 430 DEALLOCATE( & 431 & igrdi, & 432 & igrdj, & 433 & zglam, & 434 & zgphi, & 435 & zmask, & 436 & zintt, & 437 & zints & 517 DEALLOCATE( & 518 & igrdi, & 519 & igrdj, & 520 & zglam, & 521 & zgphi, & 522 & zmask, & 523 & zint, & 524 & zgdept, & 525 & zgdepw & 438 526 & ) 527 528 IF ( prodatqc%lclim ) DEALLOCATE( zclim ) 529 439 530 ! At the end of the day also get interpolated means 440 IF ( idayend == 0 ) THEN 441 DEALLOCATE( & 442 & zinmt, & 443 & zinms & 444 & ) 445 ENDIF 446 447 prodatqc%nprofup = prodatqc%nprofup + ipro 448 449 END SUBROUTINE obs_pro_opt 450 451 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 & psshn, psshmask, k2dint ) 531 IF ( ld_dailyav .AND. idayend == 0 ) THEN 532 DEALLOCATE( zinm ) 533 ENDIF 534 535 IF ( kvar == prodatqc%nvar ) THEN 536 prodatqc%nprofup = prodatqc%nprofup + ipro 537 ENDIF 538 539 END SUBROUTINE obs_prof_opt 540 541 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 542 & kit000, kdaystp, psurf, pclim, psurfmask, & 543 & k2dint, ldnightav, plamscl, pphiscl, & 544 & lindegrees, kmeanstp ) 545 453 546 !!----------------------------------------------------------------------- 454 547 !! 455 !! *** ROUTINE obs_s la_opt ***456 !! 457 !! ** Purpose : Compute the model counterpart of s ea level anomaly548 !! *** ROUTINE obs_surf_opt *** 549 !! 550 !! ** Purpose : Compute the model counterpart of surface 458 551 !! data by interpolating from the model grid to the 459 552 !! observation point. … … 462 555 !! the model values at the corners of the surrounding grid box. 463 556 !! 464 !! The n ow model SSHis first computed at the obs (lon, lat) point.557 !! The new model value is first computed at the obs (lon, lat) point. 465 558 !! 466 559 !! Several horizontal interpolation schemes are available: … … 470 563 !! - bilinear (quadrilateral grid) (k2dint = 3) 471 564 !! - polynomial (quadrilateral grid) (k2dint = 4) 472 !! 473 !! The sea level anomaly at the observation points is then computed 474 !! by removing a mean dynamic topography (defined at the obs. point). 565 !! 566 !! Two horizontal averaging schemes are also available: 567 !! - weighted radial footprint (k2dint = 5) 568 !! - weighted rectangular footprint (k2dint = 6) 569 !! 475 570 !! 476 571 !! ** Action : … … 478 573 !! History : 479 574 !! ! 07-03 (A. Weaver) 575 !! ! 15-02 (M. Martin) Combined routine for surface types 576 !! ! 17-03 (M. Martin) Added horizontal averaging options 480 577 !!----------------------------------------------------------------------- 481 578 482 579 !! * Modules used 483 580 USE obs_surf_def ! Definition of storage space for surface observations … … 486 583 487 584 !! * Arguments 488 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening 489 INTEGER, INTENT(IN) :: kt ! Time step 490 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 585 TYPE(obs_surf), INTENT(INOUT) :: & 586 & surfdataqc ! Subset of surface data passing QC 587 INTEGER, INTENT(IN) :: kt ! Time step 588 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 491 589 INTEGER, INTENT(IN) :: kpj 492 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 493 ! (kit000-1 = restart time) 494 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 495 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 496 & psshn, & ! Model SSH field 497 & psshmask ! Land-sea mask 498 590 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 591 ! (kit000-1 = restart time) 592 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 593 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 594 INTEGER, INTENT(IN), OPTIONAL :: & 595 kmeanstp ! Number of time steps for the time meaning 596 ! Averaging is triggered if present and greater than one 597 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 598 & psurf, & ! Model surface field 599 & pclim, & ! Climatological surface field 600 & psurfmask ! Land-sea mask 601 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 602 REAL(KIND=wp), INTENT(IN) :: & 603 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 604 & pphiscl ! This is the full width (rather than half-width) 605 LOGICAL, INTENT(IN) :: & 606 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 607 499 608 !! * Local declarations 500 609 INTEGER :: ji … … 502 611 INTEGER :: jobs 503 612 INTEGER :: inrc 504 INTEGER :: is la613 INTEGER :: isurf 505 614 INTEGER :: iobs 506 REAL(KIND=wp) :: zlam 507 REAL(KIND=wp) :: zphi 508 REAL(KIND=wp) :: zext(1), zobsmask(1) 509 REAL(kind=wp), DIMENSION(2,2,1) :: & 510 & zweig 615 INTEGER :: imaxifp, imaxjfp 616 INTEGER :: imodi, imodj 617 INTEGER :: idayend 618 INTEGER :: imeanend 619 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 620 & igrdi, & 621 & igrdj, & 622 & igrdip1, & 623 & igrdjp1 624 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 625 & icount_night, & 626 & imask_night 627 REAL(wp) :: zlam 628 REAL(wp) :: zphi 629 REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 630 REAL(wp) :: zdaystp 631 REAL(wp) :: zmeanstp 511 632 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 512 & zmask, & 513 & zsshl, & 514 & zglam, & 515 & zgphi 516 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 517 & igrdi, & 518 & igrdj 519 633 & zweig, & 634 & zmask, & 635 & zsurf, & 636 & zsurfm, & 637 & zsurftmp, & 638 & zclim, & 639 & zglam, & 640 & zgphi, & 641 & zglamf, & 642 & zgphif 643 644 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 645 & zintmp, & 646 & zouttmp, & 647 & zmeanday ! to compute model sst in region of 24h daylight (pole) 648 649 LOGICAL :: l_timemean 650 520 651 !------------------------------------------------------------------------ 521 652 ! Local initialization 522 653 !------------------------------------------------------------------------ 523 ! ...Record and data counters654 ! Record and data counters 524 655 inrc = kt - kit000 + 2 525 isla = sladatqc%nsstp(inrc) 656 isurf = surfdataqc%nsstp(inrc) 657 658 l_timemean = .FALSE. 659 IF ( PRESENT( kmeanstp ) ) THEN 660 IF ( kmeanstp > 1 ) l_timemean = .TRUE. 661 ENDIF 662 663 ! Work out the maximum footprint size for the 664 ! interpolation/averaging in model grid-points - has to be even. 665 666 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 667 668 669 IF ( l_timemean ) THEN 670 ! Initialize time mean for first timestep 671 imeanend = MOD( kt - kit000 + 1, kmeanstp ) 672 IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend 673 674 ! Added kt == 0 test to catch restart case 675 IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN 676 IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 677 DO jj = 1, jpj 678 DO ji = 1, jpi 679 surfdataqc%vdmean(ji,jj) = 0.0 680 END DO 681 END DO 682 ENDIF 683 684 ! On each time-step, increment the field for computing time mean 685 IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt 686 DO jj = 1, jpj 687 DO ji = 1, jpi 688 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 689 & + psurf(ji,jj) 690 END DO 691 END DO 692 693 ! Compute the time mean at the end of time period 694 IF ( imeanend == 0 ) THEN 695 zmeanstp = 1.0 / REAL( kmeanstp ) 696 IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp 697 DO jj = 1, jpj 698 DO ji = 1, jpi 699 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 700 & * zmeanstp 701 END DO 702 END DO 703 ENDIF 704 ENDIF !l_timemean 705 706 707 IF ( ldnightav ) THEN 708 709 ! Initialize array for night mean 710 IF ( kt == 0 ) THEN 711 ALLOCATE ( icount_night(kpi,kpj) ) 712 ALLOCATE ( imask_night(kpi,kpj) ) 713 ALLOCATE ( zintmp(kpi,kpj) ) 714 ALLOCATE ( zouttmp(kpi,kpj) ) 715 ALLOCATE ( zmeanday(kpi,kpj) ) 716 nday_qsr = -1 ! initialisation flag for nbc_dcy 717 ENDIF 718 719 ! Night-time means are calculated for night-time values over timesteps: 720 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 721 idayend = MOD( kt - kit000 + 1, kdaystp ) 722 723 ! Initialize night-time mean for first timestep of the day 724 IF ( idayend == 1 .OR. kt == 0 ) THEN 725 DO jj = 1, jpj 726 DO ji = 1, jpi 727 surfdataqc%vdmean(ji,jj) = 0.0 728 zmeanday(ji,jj) = 0.0 729 icount_night(ji,jj) = 0 730 END DO 731 END DO 732 ENDIF 733 734 zintmp(:,:) = 0.0 735 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 736 imask_night(:,:) = INT( zouttmp(:,:) ) 737 738 DO jj = 1, jpj 739 DO ji = 1, jpi 740 ! Increment the temperature field for computing night mean and counter 741 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 742 & + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 743 zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj) 744 icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 745 END DO 746 END DO 747 748 ! Compute the night-time mean at the end of the day 749 zdaystp = 1.0 / REAL( kdaystp ) 750 IF ( idayend == 0 ) THEN 751 IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 ! Test if "no night" point 755 IF ( icount_night(ji,jj) > 0 ) THEN 756 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 757 & / REAL( icount_night(ji,jj) ) 758 ELSE 759 !At locations where there is no night (e.g. poles), 760 ! calculate daily mean instead of night-time mean. 761 surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 762 ENDIF 763 END DO 764 END DO 765 ENDIF 766 767 ENDIF 526 768 527 769 ! Get the data for interpolation 528 770 529 771 ALLOCATE( & 530 & igrdi(2,2,isla), & 531 & igrdj(2,2,isla), & 532 & zglam(2,2,isla), & 533 & zgphi(2,2,isla), & 534 & zmask(2,2,isla), & 535 & zsshl(2,2,isla) & 772 & zweig(imaxifp,imaxjfp,1), & 773 & igrdi(imaxifp,imaxjfp,isurf), & 774 & igrdj(imaxifp,imaxjfp,isurf), & 775 & zglam(imaxifp,imaxjfp,isurf), & 776 & zgphi(imaxifp,imaxjfp,isurf), & 777 & zmask(imaxifp,imaxjfp,isurf), & 778 & zsurf(imaxifp,imaxjfp,isurf), & 779 & zsurftmp(imaxifp,imaxjfp,isurf), & 780 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 781 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 782 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 783 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 536 784 & ) 785 786 IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) ) 787 788 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 789 iobs = jobs - surfdataqc%nsurfup 790 DO ji = 0, imaxifp 791 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 792 793 !Deal with wrap around in longitude 794 IF ( imodi < 1 ) imodi = imodi + jpiglo 795 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 796 797 DO jj = 0, imaxjfp 798 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 799 !If model values are out of the domain to the north/south then 800 !set them to be the edge of the domain 801 IF ( imodj < 1 ) imodj = 1 802 IF ( imodj > jpjglo ) imodj = jpjglo 803 804 igrdip1(ji+1,jj+1,iobs) = imodi 805 igrdjp1(ji+1,jj+1,iobs) = imodj 806 807 IF ( ji >= 1 .AND. jj >= 1 ) THEN 808 igrdi(ji,jj,iobs) = imodi 809 igrdj(ji,jj,iobs) = imodj 810 ENDIF 811 812 END DO 813 END DO 814 END DO 815 816 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 817 & igrdi, igrdj, glamt, zglam ) 818 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 819 & igrdi, igrdj, gphit, zgphi ) 820 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 821 & igrdi, igrdj, psurfmask, zmask ) 822 823 ! At the end of the averaging period get interpolated means 824 IF ( l_timemean ) THEN 825 IF ( imeanend == 0 ) THEN 826 ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) ) 827 IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 828 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 829 & igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm ) 830 ENDIF 831 ELSE 832 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 833 & igrdi, igrdj, psurf, zsurf ) 834 ENDIF 835 836 IF ( k2dint > 4 ) THEN 837 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 838 & igrdip1, igrdjp1, glamf, zglamf ) 839 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 840 & igrdip1, igrdjp1, gphif, zgphif ) 841 ENDIF 537 842 538 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 539 iobs = jobs - sladatqc%nsurfup 540 igrdi(1,1,iobs) = sladatqc%mi(jobs)-1 541 igrdj(1,1,iobs) = sladatqc%mj(jobs)-1 542 igrdi(1,2,iobs) = sladatqc%mi(jobs)-1 543 igrdj(1,2,iobs) = sladatqc%mj(jobs) 544 igrdi(2,1,iobs) = sladatqc%mi(jobs) 545 igrdj(2,1,iobs) = sladatqc%mj(jobs)-1 546 igrdi(2,2,iobs) = sladatqc%mi(jobs) 547 igrdj(2,2,iobs) = sladatqc%mj(jobs) 548 END DO 549 550 CALL obs_int_comm_2d( 2, 2, isla, & 551 & igrdi, igrdj, glamt, zglam ) 552 CALL obs_int_comm_2d( 2, 2, isla, & 553 & igrdi, igrdj, gphit, zgphi ) 554 CALL obs_int_comm_2d( 2, 2, isla, & 555 & igrdi, igrdj, psshmask, zmask ) 556 CALL obs_int_comm_2d( 2, 2, isla, & 557 & igrdi, igrdj, psshn, zsshl ) 843 IF ( surfdataqc%lclim ) THEN 844 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 845 & igrdi, igrdj, pclim, zclim ) 846 ENDIF 847 848 ! At the end of the day get interpolated means 849 IF ( idayend == 0 .AND. ldnightav ) THEN 850 851 ALLOCATE( & 852 & zsurfm(imaxifp,imaxjfp,isurf) & 853 & ) 854 855 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 856 & surfdataqc%vdmean(:,:), zsurfm ) 857 858 ENDIF 558 859 559 860 ! Loop over observations 560 561 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 562 563 iobs = jobs - sladatqc%nsurfup 564 565 IF ( kt /= sladatqc%mstp(jobs) ) THEN 566 861 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 862 863 iobs = jobs - surfdataqc%nsurfup 864 865 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 866 567 867 IF(lwp) THEN 568 868 WRITE(numout,*) … … 574 874 WRITE(numout,*) ' Record = ', jobs, & 575 875 & ' kt = ', kt, & 576 & ' mstp = ', s ladatqc%mstp(jobs), &577 & ' ntyp = ', s ladatqc%ntyp(jobs)876 & ' mstp = ', surfdataqc%mstp(jobs), & 877 & ' ntyp = ', surfdataqc%ntyp(jobs) 578 878 ENDIF 579 CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 580 581 ENDIF 582 583 zlam = sladatqc%rlam(jobs) 584 zphi = sladatqc%rphi(jobs) 585 586 ! Get weights to interpolate the model SSH to the observation point 587 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 588 & zglam(:,:,iobs), zgphi(:,:,iobs), & 589 & zmask(:,:,iobs), zweig, zobsmask ) 590 591 592 ! Interpolate the model SSH to the observation point 593 CALL obs_int_h2d( 1, 1, & 594 & zweig, zsshl(:,:,iobs), zext ) 595 596 sladatqc%rext(jobs,1) = zext(1) 597 ! ... Remove the MDT at the observation point 598 sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 879 CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 880 881 ENDIF 882 883 zlam = surfdataqc%rlam(jobs) 884 zphi = surfdataqc%rphi(jobs) 885 886 IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN 887 ! Night-time or N=kmeanstp timestep averaged data 888 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 889 ELSE 890 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 891 ENDIF 892 893 IF ( ( .NOT. l_timemean ) .OR. & 894 & ( l_timemean .AND. imeanend == 0) ) THEN 895 IF ( k2dint <= 4 ) THEN 896 897 ! Get weights to interpolate the model value to the observation point 898 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 899 & zglam(:,:,iobs), zgphi(:,:,iobs), & 900 & zmask(:,:,iobs), zweig, zobsmask ) 901 902 ! Interpolate the model value to the observation point 903 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 904 905 IF ( surfdataqc%lclim ) THEN 906 CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 907 ENDIF 908 909 910 ELSE 911 912 ! Get weights to average the model SLA to the observation footprint 913 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 914 & zglam(:,:,iobs), zgphi(:,:,iobs), & 915 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 916 & zmask(:,:,iobs), plamscl, pphiscl, & 917 & lindegrees, zweig, zobsmask ) 918 919 ! Average the model SST to the observation footprint 920 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 921 & zweig, zsurftmp(:,:,iobs), zext ) 922 923 IF ( surfdataqc%lclim ) THEN 924 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 925 & zweig, zclim(:,:,iobs), zclm ) 926 ENDIF 927 928 ENDIF 929 930 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 931 ! ... Remove the MDT from the SSH at the observation point to get the SLA 932 surfdataqc%rext(jobs,1) = zext(1) 933 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 934 ELSE 935 surfdataqc%rmod(jobs,1) = zext(1) 936 ENDIF 937 938 IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 939 940 IF ( zext(1) == obfillflt ) THEN 941 ! If the observation value is a fill value, set QC flag to bad 942 surfdataqc%nqc(jobs) = 4 943 ENDIF 944 ENDIF 599 945 600 946 END DO … … 602 948 ! Deallocate the data for interpolation 603 949 DEALLOCATE( & 950 & zweig, & 604 951 & igrdi, & 605 952 & igrdj, & … … 607 954 & zgphi, & 608 955 & zmask, & 609 & zsshl & 956 & zsurf, & 957 & zsurftmp, & 958 & zglamf, & 959 & zgphif, & 960 & igrdip1,& 961 & igrdjp1 & 610 962 & ) 611 963 612 sladatqc%nsurfup = sladatqc%nsurfup + isla 613 614 END SUBROUTINE obs_sla_opt 615 616 SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 617 & psstn, psstmask, k2dint, ld_nightav ) 618 !!----------------------------------------------------------------------- 619 !! 620 !! *** ROUTINE obs_sst_opt *** 621 !! 622 !! ** Purpose : Compute the model counterpart of surface temperature 623 !! data by interpolating from the model grid to the 624 !! observation point. 625 !! 626 !! ** Method : Linearly interpolate to each observation point using 627 !! the model values at the corners of the surrounding grid box. 628