 Timestamp:
 20190701T12:44:06+02:00 (15 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4990 r11202 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 … … 52 43 & dia_obs_dealloc ! Deallocate dia_obs data 53 44 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 45 !! * Module variables 63 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles 64 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles 65 LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set 66 LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set 67 LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles 68 LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies 69 LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files 70 LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files 71 LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature 72 LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature 73 LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data 74 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files 75 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration 76 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations 77 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data 78 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data 79 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data 80 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data 81 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files 82 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height 83 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity 84 LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations 85 LOGICAL, PUBLIC :: ln_nea !: Remove observations near land 86 LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias 87 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.HHMMSS 91 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS 92 93 INTEGER, PUBLIC :: n1dint !: Vertical interpolation method 94 INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method 95 46 LOGICAL, PUBLIC :: & 47 & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. 48 LOGICAL :: ln_diaobs !: Logical switch for the obs operator 49 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 50 LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 51 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 52 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 54 LOGICAL :: ln_sic_fp_indegs !: T=> seaice obs footprint size specified in degrees, F=> in metres 55 56 REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 57 REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 58 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint 59 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint 60 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint 61 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint 62 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint 63 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint 64 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of seaice observation footprint 65 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of seaice observation footprint 66 67 INTEGER :: nn_1dint !: Vertical interpolation method 68 INTEGER :: nn_2dint_default !: Default horizontal interpolation method 69 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (1 = default) 70 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (1 = default) 71 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (1 = default) 72 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method (1 = default) 73 96 74 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? 75 & nn_profdavtypes !: Profile data types representing a daily average 76 INTEGER :: nproftypes !: Number of profile obs types 77 INTEGER :: nsurftypes !: Number of surface obs types 78 INTEGER, DIMENSION(:), ALLOCATABLE :: & 79 & nvarsprof, & !: Number of profile variables 80 & nvarssurf !: Number of surface variables 81 INTEGER, DIMENSION(:), ALLOCATABLE :: & 82 & nextrprof, & !: Number of profile extra variables 83 & nextrsurf !: Number of surface extra variables 84 INTEGER, DIMENSION(:), ALLOCATABLE :: & 85 & n2dintsurf !: Interpolation option for surface variables 86 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 87 & ravglamscl, & !: E/W diameter of averaging footprint for surface variables 88 & ravgphiscl !: N/S diameter of averaging footprint for surface variables 107 89 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 90 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 91 & llnightav !: Logical for calculating nighttime averages 92 93 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 94 & surfdata, & !: Initial surface data 95 & surfdataqc !: Surface data after quality control 96 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 97 & profdata, & !: Initial profile data 98 & profdataqc !: Profile data after quality control 99 100 CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 101 & cobstypesprof, & !: Profile obs types 102 & cobstypessurf !: Surface obs types 113 103 114 104 !! … … 118 108 !! 119 109 110 !! * Substitutions 111 # include "domzgr_substitute.h90" 120 112 CONTAINS 121 113 … … 135 127 !! ! 0610 (A. Weaver) Cleaning and add controls 136 128 !! ! 0703 (K. Mogensen) General handling of profiles 129 !! ! 1408 (J.While) Incorporated SST bias correction 130 !! ! 1502 (M. Martin) Simplification of namelist and code 137 131 !! 138 132 … … 140 134 141 135 !! * 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 136 INTEGER, PARAMETER :: & 137 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 138 INTEGER, DIMENSION(:), ALLOCATABLE :: & 139 & ifilesprof, & ! Number of profile files 140 & ifilessurf ! Number of surface files 141 INTEGER :: ios ! Local integer output status for namelist read 142 INTEGER :: jtype ! Counter for obs types 143 INTEGER :: jvar ! Counter for variables 144 INTEGER :: jfile ! Counter for files 145 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 146 INTEGER :: n2dint_type ! Local version of nn_2dint* 147 148 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 149 & cn_profbfiles, & ! T/S profile input filenames 150 & cn_sstfbfiles, & ! Sea surface temperature input filenames 151 & cn_slafbfiles, & ! Sea level anomaly input filenames 152 & cn_sicfbfiles, & ! Seaice concentration input filenames 153 & cn_velfbfiles, & ! Velocity profile input filenames 154 & cn_sssfbfiles, & ! Sea surface salinity input filenames 155 & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames 156 & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames 157 & cn_slchlnonfbfiles, & ! Surface nondiatom log10(chlorophyll) input filenames 158 & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames 159 & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 160 & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames 161 & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames 162 & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames 163 & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames 164 & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames 165 & cn_slphynonfbfiles, & ! Surface nondiatom log10(phytoplankton carbon) input filenames 166 & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames 167 & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames 168 & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames 169 & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 170 & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames 171 & cn_pno3fbfiles, & ! Profile nitrate input filenames 172 & cn_psi4fbfiles, & ! Profile silicate input filenames 173 & cn_ppo4fbfiles, & ! Profile phosphate input filenames 174 & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames 175 & cn_palkfbfiles, & ! Profile alkalinity input filenames 176 & cn_pphfbfiles, & ! Profile pH input filenames 177 & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames 178 & cn_sstbiasfiles ! SST bias input filenames 179 180 CHARACTER(LEN=128) :: & 181 & cn_altbiasfile ! Altimeter bias input filename 182 183 184 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 185 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 186 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 187 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 188 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 189 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 190 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 191 LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs 192 LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs 193 LOGICAL :: ln_slchlnon ! Logical switch for surface nondiatom log10(chlorophyll) obs 194 LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs 195 LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 196 LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs 197 LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs 198 LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs 199 LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs 200 LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs 201 LOGICAL :: ln_slphynon ! Logical switch for surface nondiatom log10(phytoplankton carbon) obs 202 LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs 203 LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs 204 LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs 205 LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs 206 LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs 207 LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs 208 LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs 209 LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs 210 LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs 211 LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs 212 LOGICAL :: ln_pph ! Logical switch for profile pH obs 213 LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs 214 LOGICAL :: ln_nea ! Logical switch to remove obs near land 215 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 216 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 217 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 218 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 219 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 220 221 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 222 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 223 224 REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 225 REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 226 227 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 228 & clproffiles, & ! Profile filenames 229 & clsurffiles ! Surface filenames 230 231 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read 232 LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 233 LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) 234 235 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 236 & zglam ! Model longitudes for profile variables 237 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 238 & zgphi ! Model latitudes for profile variables 239 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 240 & zmask ! Model land/sea mask associated with variables 241 242 243 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 244 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 245 & ln_slchltot, ln_slchldia, ln_slchlnon, & 246 & ln_slchldin, ln_slchlmic, ln_slchlnan, & 247 & ln_slchlpic, ln_schltot, & 248 & ln_slphytot, ln_slphydia, ln_slphynon, & 249 & ln_sspm, ln_sfco2, ln_spco2, & 250 & ln_plchltot, ln_pchltot, ln_pno3, & 251 & ln_psi4, ln_ppo4, ln_pdic, & 252 & ln_palk, ln_pph, ln_po2, & 253 & ln_altbias, ln_sstbias, ln_nea, & 254 & ln_grid_global, ln_grid_search_lookup, & 255 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 256 & ln_sstnight, ln_default_fp_indegs, & 257 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 258 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 259 & cn_profbfiles, cn_slafbfiles, & 260 & cn_sstfbfiles, cn_sicfbfiles, & 261 & cn_velfbfiles, cn_sssfbfiles, & 262 & cn_slchltotfbfiles, cn_slchldiafbfiles, & 263 & cn_slchlnonfbfiles, cn_slchldinfbfiles, & 264 & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & 265 & cn_slchlpicfbfiles, cn_schltotfbfiles, & 266 & cn_slphytotfbfiles, cn_slphydiafbfiles, & 267 & cn_slphynonfbfiles, cn_sspmfbfiles, & 268 & cn_sfco2fbfiles, cn_spco2fbfiles, & 269 & cn_plchltotfbfiles, cn_pchltotfbfiles, & 270 & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 271 & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & 272 & cn_po2fbfiles, & 273 & cn_sstbiasfiles, cn_altbiasfile, & 274 & cn_gridsearchfile, rn_gridsearchres, & 275 & rn_dobsini, rn_dobsend, & 276 & rn_default_avglamscl, rn_default_avgphiscl, & 277 & rn_sla_avglamscl, rn_sla_avgphiscl, & 278 & rn_sst_avglamscl, rn_sst_avgphiscl, & 279 & rn_sss_avglamscl, rn_sss_avgphiscl, & 280 & rn_sic_avglamscl, rn_sic_avgphiscl, & 281 & nn_1dint, nn_2dint_default, & 282 & nn_2dint_sla, nn_2dint_sst, & 283 & nn_2dint_sss, nn_2dint_sic, & 284 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 285 & nn_profdavtypes 205 286 206 287 ! … … 208 289 ! 209 290 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 291 ! Some namelist arrays need initialising 292 cn_profbfiles(:) = '' 293 cn_slafbfiles(:) = '' 294 cn_sstfbfiles(:) = '' 295 cn_sicfbfiles(:) = '' 296 cn_velfbfiles(:) = '' 297 cn_sssfbfiles(:) = '' 298 cn_slchltotfbfiles(:) = '' 299 cn_slchldiafbfiles(:) = '' 300 cn_slchlnonfbfiles(:) = '' 301 cn_slchldinfbfiles(:) = '' 302 cn_slchlmicfbfiles(:) = '' 303 cn_slchlnanfbfiles(:) = '' 304 cn_slchlpicfbfiles(:) = '' 305 cn_schltotfbfiles(:) = '' 306 cn_slphytotfbfiles(:) = '' 307 cn_slphydiafbfiles(:) = '' 308 cn_slphynonfbfiles(:) = '' 309 cn_sspmfbfiles(:) = '' 310 cn_sfco2fbfiles(:) = '' 311 cn_spco2fbfiles(:) = '' 312 cn_plchltotfbfiles(:) = '' 313 cn_pchltotfbfiles(:) = '' 314 cn_pno3fbfiles(:) = '' 315 cn_psi4fbfiles(:) = '' 316 cn_ppo4fbfiles(:) = '' 317 cn_pdicfbfiles(:) = '' 318 cn_palkfbfiles(:) = '' 319 cn_pphfbfiles(:) = '' 320 cn_po2fbfiles(:) = '' 321 cn_sstbiasfiles(:) = '' 322 nn_profdavtypes(:) = 1 323 324 CALL ini_date( rn_dobsini ) 325 CALL fin_date( rn_dobsend ) 326 327 ! Read namelist namobs : control observation diagnostics 328 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 240 329 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 241 330 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 242 331 243 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation332 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 244 333 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 245 334 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 246 335 IF(lwm) WRITE ( numond, namobs ) 247 336 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) 337 lk_diaobs = .FALSE. 338 #if defined key_diaobs 339 IF ( ln_diaobs ) lk_diaobs = .TRUE. 340 #endif 341 342 IF ( .NOT. lk_diaobs ) THEN 343 IF(lwp) WRITE(numout,cform_war) 344 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' 345 RETURN 253 346 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 347 322 348 IF(lwp) THEN 323 349 WRITE(numout,*) … … 325 351 WRITE(numout,*) '~~~~~~~~~~~~' 326 352 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 nighttime 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 353 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 354 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 355 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 356 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 357 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 358 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 359 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 360 WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot 361 WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia 362 WRITE(numout,*) ' Logical switch for surface nondiatom logchl obs ln_slchlnon = ', ln_slchlnon 363 WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin 364 WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic 365 WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan 366 WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic 367 WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot 368 WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot 369 WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia 370 WRITE(numout,*) ' Logical switch for surface nondiatom log(phyC) obs ln_slphynon = ', ln_slphynon 371 WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm 372 WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 373 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 374 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot 375 WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot 376 WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 377 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 378 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 379 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic 380 WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk 381 WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph 382 WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 383 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 384 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 352 385 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)) 386 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 387 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 388 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 389 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 390 WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default 391 WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla 392 WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst 393 WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss 394 WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic 395 WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl 396 WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl 397 WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs 398 WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl 399 WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl 400 WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs 401 WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl 402 WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl 403 WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs 404 WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl 405 WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl 406 WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs 407 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 408 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 409 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 410 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 411 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 412 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 413 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 414 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 415 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 416 WRITE(numout,*) ' Logical switch for nighttime SST obs ln_sstnight = ', ln_sstnight 417 ENDIF 418 ! 419 ! Set up list of observation types to be used 420 ! and the files associated with each type 421 ! 422 423 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & 424 & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & 425 & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) 426 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 427 & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 428 & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & 429 & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & 430 & ln_sfco2, ln_spco2 /) ) 431 432 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 433 IF(lwp) WRITE(numout,cform_war) 434 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 435 & ' are set to .FALSE. so turning off calls to dia_obs' 436 nwarn = nwarn + 1 437 lk_diaobs = .FALSE. 438 RETURN 439 ENDIF 440 441 IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes 442 IF ( nproftypes > 0 ) THEN 443 444 ALLOCATE( cobstypesprof(nproftypes) ) 445 ALLOCATE( ifilesprof(nproftypes) ) 446 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 447 448 jtype = 0 449 IF (ln_t3d .OR. ln_s3d) THEN 450 jtype = jtype + 1 451 cobstypesprof(jtype) = 'prof' 452 clproffiles(jtype,:) = cn_profbfiles 453 ENDIF 454 IF (ln_vel3d) THEN 455 jtype = jtype + 1 456 cobstypesprof(jtype) = 'vel' 457 clproffiles(jtype,:) = cn_velfbfiles 458 ENDIF 459 IF (ln_plchltot) THEN 460 jtype = jtype + 1 461 cobstypesprof(jtype) = 'plchltot' 462 clproffiles(jtype,:) = cn_plchltotfbfiles 463 ENDIF 464 IF (ln_pchltot) THEN 465 jtype = jtype + 1 466 cobstypesprof(jtype) = 'pchltot' 467 clproffiles(jtype,:) = cn_pchltotfbfiles 468 ENDIF 469 IF (ln_pno3) THEN 470 jtype = jtype + 1 471 cobstypesprof(jtype) = 'pno3' 472 clproffiles(jtype,:) = cn_pno3fbfiles 473 ENDIF 474 IF (ln_psi4) THEN 475 jtype = jtype + 1 476 cobstypesprof(jtype) = 'psi4' 477 clproffiles(jtype,:) = cn_psi4fbfiles 478 ENDIF 479 IF (ln_ppo4) THEN 480 jtype = jtype + 1 481 cobstypesprof(jtype) = 'ppo4' 482 clproffiles(jtype,:) = cn_ppo4fbfiles 483 ENDIF 484 IF (ln_pdic) THEN 485 jtype = jtype + 1 486 cobstypesprof(jtype) = 'pdic' 487 clproffiles(jtype,:) = cn_pdicfbfiles 488 ENDIF 489 IF (ln_palk) THEN 490 jtype = jtype + 1 491 cobstypesprof(jtype) = 'palk' 492 clproffiles(jtype,:) = cn_palkfbfiles 493 ENDIF 494 IF (ln_pph) THEN 495 jtype = jtype + 1 496 cobstypesprof(jtype) = 'pph' 497 clproffiles(jtype,:) = cn_pphfbfiles 498 ENDIF 499 IF (ln_po2) THEN 500 jtype = jtype + 1 501 cobstypesprof(jtype) = 'po2' 502 clproffiles(jtype,:) = cn_po2fbfiles 503 ENDIF 504 505 CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 506 507 ENDIF 508 509 IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes 510 IF ( nsurftypes > 0 ) THEN 511 512 ALLOCATE( cobstypessurf(nsurftypes) ) 513 ALLOCATE( ifilessurf(nsurftypes) ) 514 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 515 ALLOCATE(n2dintsurf(nsurftypes)) 516 ALLOCATE(ravglamscl(nsurftypes)) 517 ALLOCATE(ravgphiscl(nsurftypes)) 518 ALLOCATE(lfpindegs(nsurftypes)) 519 ALLOCATE(llnightav(nsurftypes)) 520 521 jtype = 0 522 IF (ln_sla) THEN 523 jtype = jtype + 1 524 cobstypessurf(jtype) = 'sla' 525 clsurffiles(jtype,:) = cn_slafbfiles 526 ENDIF 527 IF (ln_sst) THEN 528 jtype = jtype + 1 529 cobstypessurf(jtype) = 'sst' 530 clsurffiles(jtype,:) = cn_sstfbfiles 531 ENDIF 532 IF (ln_sic) THEN 533 jtype = jtype + 1 534 cobstypessurf(jtype) = 'sic' 535 clsurffiles(jtype,:) = cn_sicfbfiles 536 ENDIF 537 IF (ln_sss) THEN 538 jtype = jtype + 1 539 cobstypessurf(jtype) = 'sss' 540 clsurffiles(jtype,:) = cn_sssfbfiles 541 ENDIF 542 IF (ln_slchltot) THEN 543 jtype = jtype + 1 544 cobstypessurf(jtype) = 'slchltot' 545 clsurffiles(jtype,:) = cn_slchltotfbfiles 546 ENDIF 547 IF (ln_slchldia) THEN 548 jtype = jtype + 1 549 cobstypessurf(jtype) = 'slchldia' 550 clsurffiles(jtype,:) = cn_slchldiafbfiles 551 ENDIF 552 IF (ln_slchlnon) THEN 553 jtype = jtype + 1 554 cobstypessurf(jtype) = 'slchlnon' 555 clsurffiles(jtype,:) = cn_slchlnonfbfiles 556 ENDIF 557 IF (ln_slchldin) THEN 558 jtype = jtype + 1 559 cobstypessurf(jtype) = 'slchldin' 560 clsurffiles(jtype,:) = cn_slchldinfbfiles 561 ENDIF 562 IF (ln_slchlmic) THEN 563 jtype = jtype + 1 564 cobstypessurf(jtype) = 'slchlmic' 565 clsurffiles(jtype,:) = cn_slchlmicfbfiles 566 ENDIF 567 IF (ln_slchlnan) THEN 568 jtype = jtype + 1 569 cobstypessurf(jtype) = 'slchlnan' 570 clsurffiles(jtype,:) = cn_slchlnanfbfiles 571 ENDIF 572 IF (ln_slchlpic) THEN 573 jtype = jtype + 1 574 cobstypessurf(jtype) = 'slchlpic' 575 clsurffiles(jtype,:) = cn_slchlpicfbfiles 576 ENDIF 577 IF (ln_schltot) THEN 578 jtype = jtype + 1 579 cobstypessurf(jtype) = 'schltot' 580 clsurffiles(jtype,:) = cn_schltotfbfiles 581 ENDIF 582 IF (ln_slphytot) THEN 583 jtype = jtype + 1 584 cobstypessurf(jtype) = 'slphytot' 585 clsurffiles(jtype,:) = cn_slphytotfbfiles 586 ENDIF 587 IF (ln_slphydia) THEN 588 jtype = jtype + 1 589 cobstypessurf(jtype) = 'slphydia' 590 clsurffiles(jtype,:) = cn_slphydiafbfiles 591 ENDIF 592 IF (ln_slphynon) THEN 593 jtype = jtype + 1 594 cobstypessurf(jtype) = 'slphynon' 595 clsurffiles(jtype,:) = cn_slphynonfbfiles 596 ENDIF 597 IF (ln_sspm) THEN 598 jtype = jtype + 1 599 cobstypessurf(jtype) = 'sspm' 600 clsurffiles(jtype,:) = cn_sspmfbfiles 601 ENDIF 602 IF (ln_sfco2) THEN 603 jtype = jtype + 1 604 cobstypessurf(jtype) = 'sfco2' 605 clsurffiles(jtype,:) = cn_sfco2fbfiles 606 ENDIF 607 IF (ln_spco2) THEN 608 jtype = jtype + 1 609 cobstypessurf(jtype) = 'spco2' 610 clsurffiles(jtype,:) = cn_spco2fbfiles 611 ENDIF 612 613 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 614 615 DO jtype = 1, nsurftypes 616 617 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 618 IF ( nn_2dint_sla == 1 ) THEN 619 n2dint_type = nn_2dint_default 371 620 ELSE 372 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 373 TRIM(profbfiles(ji)) 621 n2dint_type = nn_2dint_sla 374 622 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)) 623 ztype_avglamscl = rn_sla_avglamscl 624 ztype_avgphiscl = rn_sla_avgphiscl 625 ltype_fp_indegs = ln_sla_fp_indegs 626 ltype_night = .FALSE. 627 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 628 IF ( nn_2dint_sst == 1 ) THEN 629 n2dint_type = nn_2dint_default 441 630 ELSE 442 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 443 TRIM(velfbfiles(ji)) 631 n2dint_type = nn_2dint_sst 444 632 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 633 ztype_avglamscl = rn_sst_avglamscl 634 ztype_avgphiscl = rn_sst_avgphiscl 635 ltype_fp_indegs = ln_sst_fp_indegs 636 ltype_night = ln_sstnight 637 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 638 IF ( nn_2dint_sic == 1 ) THEN 639 n2dint_type = nn_2dint_default 640 ELSE 641 n2dint_type = nn_2dint_sic 642 ENDIF 643 ztype_avglamscl = rn_sic_avglamscl 644 ztype_avgphiscl = rn_sic_avgphiscl 645 ltype_fp_indegs = ln_sic_fp_indegs 646 ltype_night = .FALSE. 647 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 648 IF ( nn_2dint_sss == 1 ) THEN 649 n2dint_type = nn_2dint_default 650 ELSE 651 n2dint_type = nn_2dint_sss 652 ENDIF 653 ztype_avglamscl = rn_sss_avglamscl 654 ztype_avgphiscl = rn_sss_avgphiscl 655 ltype_fp_indegs = ln_sss_fp_indegs 656 ltype_night = .FALSE. 657 ELSE 658 n2dint_type = nn_2dint_default 659 ztype_avglamscl = rn_default_avglamscl 660 ztype_avgphiscl = rn_default_avgphiscl 661 ltype_fp_indegs = ln_default_fp_indegs 662 ltype_night = .FALSE. 663 ENDIF 664 665 CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 666 & nn_2dint_default, n2dint_type, & 667 & ztype_avglamscl, ztype_avgphiscl, & 668 & ltype_fp_indegs, ltype_night, & 669 & n2dintsurf, ravglamscl, ravgphiscl, & 670 & lfpindegs, llnightav ) 671 672 END DO 458 673 459 674 ENDIF 460 675 676 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 677 678 679 ! 680 ! Obs operator parameter checking and initialisations 681 ! 682 461 683 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 462 684 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 464 686 ENDIF 465 687 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 688 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 485 689 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 486 690 & ' is not available') 487 691 ENDIF 488 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 489 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 692 693 IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN 694 CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & 490 695 & ' is not available') 491 696 ENDIF 492 697 698 CALL obs_typ_init 699 700 CALL mppmap_init 701 702 CALL obs_grid_setup( ) 703 493 704 ! 494 705 ! Depending on switches read the various observation types 495 706 ! 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 707 708 IF ( nproftypes > 0 ) THEN 709 710 ALLOCATE(profdata(nproftypes)) 711 ALLOCATE(profdataqc(nproftypes)) 712 ALLOCATE(nvarsprof(nproftypes)) 713 ALLOCATE(nextrprof(nproftypes)) 714 715 DO jtype = 1, nproftypes 716 717 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 718 nvarsprof(jtype) = 2 719 nextrprof(jtype) = 1 720 ALLOCATE(llvar(nvarsprof(jtype))) 721 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 722 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 723 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 724 llvar(1) = ln_t3d 725 llvar(2) = ln_s3d 726 zglam(:,:,1) = glamt(:,:) 727 zglam(:,:,2) = glamt(:,:) 728 zgphi(:,:,1) = gphit(:,:) 729 zgphi(:,:,2) = gphit(:,:) 730 zmask(:,:,:,1) = tmask(:,:,:) 731 zmask(:,:,:,2) = tmask(:,:,:) 732 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 733 nvarsprof(jtype) = 2 734 nextrprof(jtype) = 2 735 ALLOCATE(llvar(nvarsprof(jtype))) 736 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 737 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 738 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 739 llvar(1) = ln_vel3d 740 llvar(2) = ln_vel3d 741 zglam(:,:,1) = glamu(:,:) 742 zglam(:,:,2) = glamv(:,:) 743 zgphi(:,:,1) = gphiu(:,:) 744 zgphi(:,:,2) = gphiv(:,:) 745 zmask(:,:,:,1) = umask(:,:,:) 746 zmask(:,:,:,2) = vmask(:,:,:) 747 ELSE 748 nvarsprof(jtype) = 1 749 nextrprof(jtype) = 0 750 ALLOCATE(llvar(nvarsprof(jtype))) 751 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 752 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 753 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 754 llvar(1) = .TRUE. 755 zglam(:,:,1) = glamt(:,:) 756 zgphi(:,:,1) = gphit(:,:) 757 zmask(:,:,:,1) = tmask(:,:,:) 758 ENDIF 759 760 !Read in profile or profile obs types 761 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 762 & clproffiles(jtype,1:ifilesprof(jtype)), & 763 & nvarsprof(jtype), nextrprof(jtype), nitendnit000+2, & 764 & rn_dobsini, rn_dobsend, llvar, & 765 & ln_ignmis, ln_s_at_t, .FALSE., & 766 & kdailyavtypes = nn_profdavtypes ) 767 768 DO jvar = 1, nvarsprof(jtype) 769 CALL obs_prof_staend( profdata(jtype), jvar ) 770 END DO 771 772 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 773 & llvar, & 774 & jpi, jpj, jpk, & 775 & zmask, zglam, zgphi, & 776 & ln_nea, ln_bound_reject, & 777 & kdailyavtypes = nn_profdavtypes ) 524 778 525 ld_enact(jprofset) = .TRUE. 526 527 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 528 & jnumenact, enactfiles(1:jnumenact), & 529 & nprofvars, nprofextr, & 530 & nitendnit000+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 539 END DO 540 541 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 542 & ln_t3d, ln_s3d, ln_nea, & 543 & kdailyavtypes=endailyavtypes ) 544 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 & nitendnit000+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 ) 570 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 & nitendnit000+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 779 DEALLOCATE( llvar ) 780 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zglam ) 781 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zgphi ) 782 CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 783 784 END DO 785 786 DEALLOCATE( ifilesprof, clproffiles ) 787 788 ENDIF 789 790 IF ( nsurftypes > 0 ) THEN 791 792 ALLOCATE(surfdata(nsurftypes)) 793 ALLOCATE(surfdataqc(nsurftypes)) 794 ALLOCATE(nvarssurf(nsurftypes)) 795 ALLOCATE(nextrsurf(nsurftypes)) 796 797 DO jtype = 1, nsurftypes 798 799 nvarssurf(jtype) = 1 800 nextrsurf(jtype) = 0 801 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 802 803 !Read in surface obs types 804 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 805 & clsurffiles(jtype,1:ifilessurf(jtype)), & 806 & nvarssurf(jtype), nextrsurf(jtype), nitendnit000+2, & 807 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 808 809 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 810 811 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 812 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 813 IF ( ln_altbias ) & 814 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 815 ENDIF 816 817 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 818 jnumsstbias = 0 819 DO jfile = 1, jpmaxnfiles 820 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 821 & jnumsstbias = jnumsstbias + 1 596 822 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 ) 823 IF ( jnumsstbias == 0 ) THEN 824 CALL ctl_stop("ln_sstbias set but no bias files to read in") 605 825 ENDIF 606 607 END DO 608 609 ENDIF 826 827 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 828 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 829 830 ENDIF 831 832 END DO 833 834 DEALLOCATE( ifilessurf, clsurffiles ) 610 835 611 836 ENDIF 612 837 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, nitendnit000+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, nitendnit000+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, nitendnit000+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 & nitendnit000+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, nitendnit000+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, nitendnit000+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, nitendnit000+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 & nitendnit000+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 & nitendnit000+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 & nitendnit000+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 & nitendnit000+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 & nitendnit000+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 838 END SUBROUTINE dia_obs_init 968 839 … … 974 845 !! 975 846 !! ** 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 983 !! 984 !! ** Action : 847 !! compute the model equivalent of the following data: 848 !!  Profile data, currently T/S or U/V 849 !!  Surface data, currently SST, SLA or seaice concentration. 850 !! 851 !! ** Action : 985 852 !! 986 853 !! History : … … 991 858 !! ! 0704 (G. Smith) Generalized surface operators 992 859 !! ! 0810 (M. Valdivieso) obs operator for velocity profiles 860 !! ! 1508 (M. Martin) Combined surface/profile routines. 993 861 !! 994 862 !! * 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,&863 USE phycst, ONLY : & ! Physical constants 864 #if defined key_fabm 865 & rt0, & 866 #endif 867 & rday 868 USE oce, ONLY : & ! Ocean dynamics and tracers variables 869 & tsn, & 870 & un, & 871 & vn, & 1004 872 & sshn 1005 873 #if defined key_lim3 1006 USE ice, ONLY : & ! LIMIce model variables874 USE ice, ONLY : & ! LIM3 Ice model variables 1007 875 & frld 1008 876 #endif 1009 877 #if defined key_lim2 1010 USE ice_2, ONLY : & ! LIMIce model variables878 USE ice_2, ONLY : & ! LIM2 Ice model variables 1011 879 & frld 1012 880 #endif 881 #if defined key_cice 882 USE sbc_oce, ONLY : fr_i ! ice fraction 883 #endif 884 #if defined key_top 885 USE trc, ONLY : & ! Biogeochemical state variables 886 & trn 887 #endif 888 #if defined key_hadocc 889 USE par_hadocc ! HadOCC parameters 890 USE trc, ONLY : & 891 & HADOCC_CHL, & 892 & HADOCC_FCO2, & 893 & HADOCC_PCO2, & 894 & HADOCC_FILL_FLT 895 USE had_bgc_const, ONLY: c2n_p 896 #elif defined key_medusa 897 USE par_medusa ! MEDUSA parameters 898 USE sms_medusa, ONLY: & 899 & xthetapn, & 900 & xthetapd 901 #if defined key_roam 902 USE sms_medusa, ONLY: & 903 & f2_pco2w, & 904 & f2_fco2w, & 905 & f3_pH 906 #endif 907 #elif defined key_fabm 908 USE par_fabm ! FABM parameters 909 USE fabm, ONLY: & 910 & fabm_get_interior_diagnostic_data 911 #endif 912 #if defined key_spm 913 USE par_spm, ONLY: & ! Sediment parameters 914 & jp_spm 915 #endif 916 1013 917 IMPLICIT NONE 1014 918 1015 919 !! * Arguments 1016 INTEGER, INTENT(IN) :: kstp 920 INTEGER, INTENT(IN) :: kstp ! Current timestep 1017 921 !! * 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) 922 INTEGER :: idaystp ! Number of timesteps per day 923 INTEGER :: jtype ! Data loop variable 924 INTEGER :: jvar ! Variable number 925 INTEGER :: ji, jj, jk ! Loop counters 926 REAL(wp) :: tiny ! small number 927 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 928 & zprofvar ! Model values for variables in a prof ob 929 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 930 & zprofmask ! Mask associated with zprofvar 931 REAL(wp), POINTER, DIMENSION(:,:) :: & 932 & zsurfvar, & ! Model values equivalent to surface ob. 933 & zsurfmask ! Mask associated with surface variable 934 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 935 & zglam, & ! Model longitudes for prof variables 936 & zgphi ! Model latitudes for prof variables 937 LOGICAL :: llog10 ! Perform log10 transform of variable 938 #if defined key_fabm 939 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 940 & pco2_3d ! 3D pCO2 from FABM 1032 941 #endif 1033 942 … … 1036 945 WRITE(numout,*) 'dia_obs : Call the observation operators', kstp 1037 946 WRITE(numout,*) '~~~~~~~' 947 CALL FLUSH(numout) 1038 948 ENDIF 1039 949 … … 1041 951 1042 952 ! 1043 ! No LIM => frld == 0.0_wp 1044 ! 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 ) 1061 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 ) 953 ! Call the profile and surface observation operators 954 ! 955 956 IF ( nproftypes > 0 ) THEN 957 958 DO jtype = 1, nproftypes 959 960 ! Allocate local work arrays 961 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 962 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 963 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 964 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 965 966 ! Defaults which might change 967 DO jvar = 1, profdataqc(jtype)%nvar 968 zprofmask(:,:,:,jvar) = tmask(:,:,:) 969 zglam(:,:,jvar) = glamt(:,:) 970 zgphi(:,:,jvar) = gphit(:,:) 971 END DO 972 973 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 974 975 CASE('prof') 976 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 977 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 978 979 CASE('vel') 980 zprofvar(:,:,:,1) = un(:,:,:) 981 zprofvar(:,:,:,2) = vn(:,:,:) 982 zprofmask(:,:,:,1) = umask(:,:,:) 983 zprofmask(:,:,:,2) = vmask(:,:,:) 984 zglam(:,:,1) = glamu(:,:) 985 zglam(:,:,2) = glamv(:,:) 986 zgphi(:,:,1) = gphiu(:,:) 987 zgphi(:,:,2) = gphiv(:,:) 988 989 CASE('plchltot') 990 #if defined key_hadocc 991 ! Chlorophyll from HadOCC 992 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 993 #elif defined key_medusa 994 ! Add nondiatom and diatom chlorophyll from MEDUSA 995 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 996 #elif defined key_fabm 997 ! Add all chlorophyll groups from ERSEM 998 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 999 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1000 #else 1001 CALL ctl_stop( ' Trying to run plchltot observation operator', & 1002 & ' but no biogeochemical model appears to have been defined' ) 1003 #endif 1004 ! Take the log10 where we can, otherwise exclude 1005 tiny = 1.0e20 1006 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 1007 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 1008 ELSEWHERE 1009 zprofvar(:,:,:,:) = obfillflt 1010 zprofmask(:,:,:,:) = 0 1011 END WHERE 1012 ! Mask out model below any excluded values, 1013 ! to avoid interpolation issues 1014 DO jvar = 1, profdataqc(jtype)%nvar 1015 DO jj = 1, jpj 1016 DO ji = 1, jpi 1017 depth_loop: DO jk = 1, jpk 1018 IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 1019 zprofmask(ji,jj,jk:jpk,jvar) = 0 1020 EXIT depth_loop 1021 ENDIF 1022 END DO depth_loop 1023 END DO 1024 END DO 1025 END DO 1026 1027 CASE('pchltot') 1028 #if defined key_hadocc 1029 ! Chlorophyll from HadOCC 1030 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1031 #elif defined key_medusa 1032 ! Add nondiatom and diatom chlorophyll from MEDUSA 1033 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1034 #elif defined key_fabm 1035 ! Add all chlorophyll groups from ERSEM 1036 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1037 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1038 #else 1039 CALL ctl_stop( ' Trying to run pchltot observation operator', & 1040 & ' but no biogeochemical model appears to have been defined' ) 1041 #endif 1042 1043 CASE('pno3') 1044 #if defined key_hadocc 1045 ! Dissolved inorganic nitrogen from HadOCC 1046 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1047 #elif defined key_medusa 1048 ! Dissolved inorganic nitrogen from MEDUSA 1049 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1050 #elif defined key_fabm 1051 ! Nitrate from ERSEM 1052 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 1053 #else 1054 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1055 & ' but no biogeochemical model appears to have been defined' ) 1056 #endif 1057 1058 CASE('psi4') 1059 #if defined key_hadocc 1060 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1061 & ' but HadOCC does not simulate silicate' ) 1062 #elif defined key_medusa 1063 ! Silicate from MEDUSA 1064 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1065 #elif defined key_fabm 1066 ! Silicate from ERSEM 1067 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 1068 #else 1069 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1070 & ' but no biogeochemical model appears to have been defined' ) 1071 #endif 1072 1073 CASE('ppo4') 1074 #if defined key_hadocc 1075 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1076 & ' but HadOCC does not simulate phosphate' ) 1077 #elif defined key_medusa 1078 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1079 & ' but MEDUSA does not simulate phosphate' ) 1080 #elif defined key_fabm 1081 ! Phosphate from ERSEM 1082 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 1083 #else 1084 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1085 & ' but no biogeochemical model appears to have been defined' ) 1086 #endif 1087 1088 CASE('pdic') 1089 #if defined key_hadocc 1090 ! Dissolved inorganic carbon from HadOCC 1091 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1092 #elif defined key_medusa 1093 ! Dissolved inorganic carbon from MEDUSA 1094 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1095 #elif defined key_fabm 1096 ! Dissolved inorganic carbon from ERSEM 1097 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 1098 #else 1099 CALL ctl_stop( ' Trying to run pdic observation operator', & 1100 & ' but no biogeochemical model appears to have been defined' ) 1101 #endif 1102 1103 CASE('palk') 1104 #if defined key_hadocc 1105 ! Alkalinity from HadOCC 1106 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1107 #elif defined key_medusa 1108 ! Alkalinity from MEDUSA 1109 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1110 #elif defined key_fabm 1111 ! Alkalinity from ERSEM 1112 zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3a) 1113 #else 1114 CALL ctl_stop( ' Trying to run palk observation operator', & 1115 & ' but no biogeochemical model appears to have been defined' ) 1116 #endif 1117 1118 CASE('pph') 1119 #if defined key_hadocc 1120 CALL ctl_stop( ' Trying to run pph observation operator', & 1121 & ' but HadOCC has no pH diagnostic defined' ) 1122 #elif defined key_medusa && defined key_roam 1123 ! pH from MEDUSA 1124 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1125 #elif defined key_fabm 1126 ! pH from ERSEM 1127 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 1128 #else 1129 CALL ctl_stop( ' Trying to run pph observation operator', & 1130 & ' but no biogeochemical model appears to have been defined' ) 1131 #endif 1132 1133 CASE('po2') 1134 #if defined key_hadocc 1135 CALL ctl_stop( ' Trying to run po2 observation operator', & 1136 & ' but HadOCC does not simulate oxygen' ) 1137 #elif defined key_medusa 1138 ! Oxygen from MEDUSA 1139 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1140 #elif defined key_fabm 1141 ! Oxygen from ERSEM 1142 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 1143 #else 1144 CALL ctl_stop( ' Trying to run po2 observation operator', & 1145 & ' but no biogeochemical model appears to have been defined' ) 1146 #endif 1147 1148 CASE DEFAULT 1149 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1150 1151 END SELECT 1152 1153 DO jvar = 1, profdataqc(jtype)%nvar 1154 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 1155 & nit000, idaystp, jvar, & 1156 & zprofvar(:,:,:,jvar), & 1157 & fsdept(:,:,:), fsdepw(:,:,:), & 1158 & zprofmask(:,:,:,jvar), & 1159 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1160 & nn_1dint, nn_2dint_default, & 1161 & kdailyavtypes = nn_profdavtypes ) 1162 END DO 1163 1164 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1165 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1166 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1167 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1168 1169 END DO 1170 1171 ENDIF 1172 1173 IF ( nsurftypes > 0 ) THEN 1174 1175 !Allocate local work arrays 1176 CALL wrk_alloc( jpi, jpj, zsurfvar ) 1177 CALL wrk_alloc( jpi, jpj, zsurfmask ) 1178 #if defined key_fabm 1179 CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 1180 #endif 1181 1182 DO jtype = 1, nsurftypes 1183 1184 !Defaults which might be changed 1185 zsurfmask(:,:) = tmask(:,:,1) 1186 llog10 = .FALSE. 1187 1188 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 1189 CASE('sst') 1190 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 1191 CASE('sla') 1192 zsurfvar(:,:) = sshn(:,:) 1193 CASE('sss') 1194 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 1195 CASE('sic') 1196 IF ( kstp == 0 ) THEN 1197 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1198 CALL ctl_warn( 'Seaice not initialised on zeroth '// & 1199 & 'timestep but some obs are valid then.' ) 1200 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1201 & ' seaice obs will be missed' 1202 ENDIF 1203 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 1204 & surfdataqc(jtype)%nsstp(1) 1205 CYCLE 1206 ELSE 1207 #if defined key_cice 1208 zsurfvar(:,:) = fr_i(:,:) 1209 #elif defined key_lim2  defined key_lim3 1210 zsurfvar(:,:) = 1._wp  frld(:,:) 1211 #else 1212 CALL ctl_stop( ' Trying to run seaice observation operator', & 1213 & ' but no seaice model appears to have been defined' ) 1214 #endif 1215 ENDIF 1216 1217 CASE('slchltot') 1218 #if defined key_hadocc 1219 ! Surface chlorophyll from HadOCC 1220 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1221 #elif defined key_medusa 1222 ! Add nondiatom and diatom surface chlorophyll from MEDUSA 1223 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1224 #elif defined key_fabm 1225 ! Add all surface chlorophyll groups from ERSEM 1226 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1227 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1228 #else 1229 CALL ctl_stop( ' Trying to run slchltot observation operator', & 1230 & ' but no biogeochemical model appears to have been defined' ) 1231 #endif 1232 llog10 = .TRUE. 1233 1234 CASE('slchldia') 1235 #if defined key_hadocc 1236 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1237 & ' but HadOCC does not explicitly simulate diatoms' ) 1238 #elif defined key_medusa 1239 ! Diatom surface chlorophyll from MEDUSA 1240 zsurfvar(:,:) = trn(:,:,1,jpchd) 1241 #elif defined key_fabm 1242 ! Diatom surface chlorophyll from ERSEM 1243 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 1244 #else 1245 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1246 & ' but no biogeochemical model appears to have been defined' ) 1247 #endif 1248 llog10 = .TRUE. 1249 1250 CASE('slchlnon') 1251 #if defined key_hadocc 1252 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1253 & ' but HadOCC does not explicitly simulate nondiatoms' ) 1254 #elif defined key_medusa 1255 ! Nondiatom surface chlorophyll from MEDUSA 1256 zsurfvar(:,:) = trn(:,:,1,jpchn) 1257 #elif defined key_fabm 1258 ! Add all nondiatom surface chlorophyll groups from ERSEM 1259 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1260 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1261 #else 1262 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1263 & ' but no biogeochemical model appears to have been defined' ) 1264 #endif 1265 llog10 = .TRUE. 1266 1267 CASE('slchldin') 1268 #if defined key_hadocc 1269 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1270 & ' but HadOCC does not explicitly simulate dinoflagellates' ) 1271 #elif defined key_medusa 1272 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1273 & ' but MEDUSA does not explicitly simulate dinoflagellates' ) 1274 #elif defined key_fabm 1275 ! Dinoflagellate surface chlorophyll from ERSEM 1276 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1277 #else 1278 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1279 & ' but no biogeochemical model appears to have been defined' ) 1280 #endif 1281 llog10 = .TRUE. 1282 1283 CASE('slchlmic') 1284 #if defined key_hadocc 1285 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1286 & ' but HadOCC does not explicitly simulate microphytoplankton' ) 1287 #elif defined key_medusa 1288 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1289 & ' but MEDUSA does not explicitly simulate microphytoplankton' ) 1290 #elif defined key_fabm 1291 ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 1292 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1293 #else 1294 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1295 & ' but no biogeochemical model appears to have been defined' ) 1296 #endif 1297 llog10 = .TRUE. 1298 1299 CASE('slchlnan') 1300 #if defined key_hadocc 1301 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1302 & ' but HadOCC does not explicitly simulate nanophytoplankton' ) 1303 #elif defined key_medusa 1304 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1305 & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 1306 #elif defined key_fabm 1307 ! Nanophytoplankton surface chlorophyll from ERSEM 1308 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 1309 #else 1310 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1311 & ' but no biogeochemical model appears to have been defined' ) 1312 #endif 1313 llog10 = .TRUE. 1314 1315 CASE('slchlpic') 1316 #if defined key_hadocc 1317 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1318 & ' but HadOCC does not explicitly simulate picophytoplankton' ) 1319 #elif defined key_medusa 1320 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1321 & ' but MEDUSA does not explicitly simulate picophytoplankton' ) 1322 #elif defined key_fabm 1323 ! Picophytoplankton surface chlorophyll from ERSEM 1324 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 1325 #else 1326 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1327 & ' but no biogeochemical model appears to have been defined' ) 1328 #endif 1329 llog10 = .TRUE. 1330 1331 CASE('schltot') 1332 #if defined key_hadocc 1333 ! Surface chlorophyll from HadOCC 1334 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1335 #elif defined key_medusa 1336 ! Add nondiatom and diatom surface chlorophyll from MEDUSA 1337 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1338 #elif defined key_fabm 1339 ! Add all surface chlorophyll groups from ERSEM 1340 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1341 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1342 #else 1343 CALL ctl_stop( ' Trying to run schltot observation operator', & 1344 & ' but no biogeochemical model appears to have been defined' ) 1345 #endif 1346 1347 CASE('slphytot') 1348 #if defined key_hadocc 1349 ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 1350 zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 1351 #elif defined key_medusa 1352 ! Add nondiatom and diatom surface phytoplankton nitrogen from MEDUSA 1353 ! multiplied by C:N ratio for each 1354 zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 1355 #elif defined key_fabm 1356 ! Add all surface phytoplankton carbon groups from ERSEM 1357 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1358 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1359 #else 1360 CALL ctl_stop( ' Trying to run slphytot observation operator', & 1361 & ' but no biogeochemical model appears to have been defined' ) 1362 #endif 1363 llog10 = .TRUE. 1364 1365 CASE('slphydia') 1366 #if defined key_hadocc 1367 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1368 & ' but HadOCC does not explicitly simulate diatoms' ) 1369 #elif defined key_medusa 1370 ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1371 zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 1372 #elif defined key_fabm 1373 ! Diatom surface phytoplankton carbon from ERSEM 1374 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 1375 #else 1376 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1377 & ' but no biogeochemical model appears to have been defined' ) 1378 #endif 1379 llog10 = .TRUE. 1380 1381 CASE('slphynon') 1382 #if defined key_hadocc 1383 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1384 & ' but HadOCC does not explicitly simulate nondiatoms' ) 1385 #elif defined key_medusa 1386 ! Nondiatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1387 zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 1388 #elif defined key_fabm 1389 ! Add all nondiatom surface phytoplankton carbon groups from ERSEM 1390 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1391 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1392 #else 1393 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1394 & ' but no biogeochemical model appears to have been defined' ) 1395 #endif 1396 llog10 = .TRUE. 1397 1398 CASE('sspm') 1399 #if defined key_spm 1400 zsurfvar(:,:) = 0.0 1401 DO jn = 1, jp_spm 1402 zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes 1403 END DO 1404 #else 1405 CALL ctl_stop( ' Trying to run sspm observation operator', & 1406 & ' but no spm model appears to have been defined' ) 1407 #endif 1408 1409 CASE('sfco2') 1410 #if defined key_hadocc 1411 zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 1412 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 1413 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 1414 zsurfvar(:,:) = obfillflt 1415 zsurfmask(:,:) = 0 1416 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 1417 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 1418 ENDIF 1419 #elif defined key_medusa && defined key_roam 1420 zsurfvar(:,:) = f2_fco2w(:,:) 1421 #elif defined key_fabm 1422 ! First, get pCO2 from FABM 1423 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1424 zsurfvar(:,:) = pco2_3d(:,:,1) 1425 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 1426 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 1427 ! and data reduction routines, DeepSea Research II, 56: 512522. 1428 ! and 1429 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a nonideal gas, 1430 ! Marine Chemistry, 2: 203215. 1431 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 1432 ! not explicitly included  atmospheric pressure is not necessarily available so this is 1433 ! the best assumption. 1434 ! Further, the (1xCO2)^2 term has been neglected. This is common practice 1435 ! (see e.g. Zeebe and WolfGladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 1436 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 1437 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 1438 zsurfvar(:,:) = zsurfvar(:,:) * EXP((1636.75 + & 1439 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0)  & 1440 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1441 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1442 & 2.0 * (57.7  0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 1443 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 1444 #else 1445 CALL ctl_stop( ' Trying to run sfco2 observation operator', & 1446 & ' but no biogeochemical model appears to have been defined' ) 1447 #endif 1448 1449 CASE('spco2') 1450 #if defined key_hadocc 1451 zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 1452 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 1453 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 1454 zsurfvar(:,:) = obfillflt 1455 zsurfmask(:,:) = 0 1456 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 1457 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 1458 ENDIF 1459 #elif defined key_medusa && defined key_roam 1460 zsurfvar(:,:) = f2_pco2w(:,:) 1461 #elif defined key_fabm 1462 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1463 zsurfvar(:,:) = pco2_3d(:,:,1) 1464 #else 1465 CALL ctl_stop( ' Trying to run spco2 observation operator', & 1466 & ' but no biogeochemical model appears to have been defined' ) 1467 #endif 1468 1469 CASE DEFAULT 1470 1471 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 1472 1473 END SELECT 1474 1475 IF ( llog10 ) THEN 1476 ! Take the log10 where we can, otherwise exclude 1477 tiny = 1.0e20 1478 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 1479 zsurfvar(:,:) = LOG10(zsurfvar(:,:)) 1480 ELSEWHERE 1481 zsurfvar(:,:) = obfillflt 1482 zsurfmask(:,:) = 0 1483 END WHERE 1066 1484 ENDIF 1485 1486 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1487 & nit000, idaystp, zsurfvar, zsurfmask, & 1488 & n2dintsurf(jtype), llnightav(jtype), & 1489 & ravglamscl(jtype), ravgphiscl(jtype), & 1490 & lfpindegs(jtype) ) 1491 1067 1492 END DO 1493 1494 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 1495 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 1496 #if defined key_fabm 1497 CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 1498 #endif 1499 1068 1500 ENDIF 1069 1501 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 1502 END SUBROUTINE dia_obs 1119 1120 SUBROUTINE dia_obs_wri 1503 1504 SUBROUTINE dia_obs_wri 1121 1505 !! 1122 1506 !! *** ROUTINE dia_obs_wri *** … … 1126 1510 !! ** Method : Call observation diagnostic output routines 1127 1511 !! 1128 !! ** Action : 1512 !! ** Action : 1129 1513 !! 1130 1514 !! History : … … 1134 1518 !! ! 0703 (K. Mogensen) General handling of profiles 1135 1519 !! ! 0809 (M. Valdivieso) Velocity component (U,V) profiles 1520 !! ! 1508 (M. Martin) Combined writing for prof and surf types 1136 1521 !! 1522 !! * Modules used 1523 USE obs_rot_vel ! Rotation of velocities 1524 1137 1525 IMPLICIT NONE 1138 1526 1139 1527 !! * 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 1528 INTEGER :: jtype ! Data set loop variable 1529 INTEGER :: jo, jvar, jk 1530 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 1531 & zu, & 1532 & zv 1533 1150 1534 ! 1151 1535 ! Depending on switches call various observation output routines 1152 1536 ! 1153 1537 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 ) 1538 IF ( nproftypes > 0 ) THEN 1539 1540 DO jtype = 1, nproftypes 1541 1542 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 1543 1544 ! For velocity data, rotate the model velocities to N/S, E/W 1545 ! using the compressed data structure. 1546 ALLOCATE( & 1547 & zu(profdataqc(jtype)%nvprot(1)), & 1548 & zv(profdataqc(jtype)%nvprot(2)) & 1549 & ) 1550 1551 CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 1552 1553 DO jo = 1, profdataqc(jtype)%nprof 1554 DO jvar = 1, 2 1555 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 1556 1557 IF ( jvar == 1 ) THEN 1558 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 1559 ELSE 1560 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 1561 ENDIF 1562 1563 END DO 1564 END DO 1565 END DO 1566 1567 DEALLOCATE( zu ) 1568 DEALLOCATE( zv ) 1569 1570 END IF 1571 1572 CALL obs_prof_decompress( profdataqc(jtype), & 1573 & profdata(jtype), .TRUE., numout ) 1574 1575 CALL obs_wri_prof( profdata(jtype) ) 1163 1576 1164 1577 END DO 1165 1578 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 1579 ENDIF 1207 1580 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)1581 IF ( nsurftypes > 0 ) THEN 1582 1583 DO jtype = 1, nsurftypes 1584 1585 CALL obs_surf_decompress( surfdataqc(jtype), & 1586 & surfdata(jtype), .TRUE., numout ) 1587 1588 CALL obs_wri_surf( surfdata(jtype) ) 1216 1589 1217 1590 END DO 1218 1591 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 1592 ENDIF 1391 1593 … … 1405 1607 !! 1406 1608 !! 1407 ! !obs_grid deallocation1609 ! obs_grid deallocation 1408 1610 CALL obs_grid_deallocate 1409 1611 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 1612 ! diaobs deallocation 1613 IF ( nproftypes > 0 ) & 1614 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 1615 1616 IF ( nsurftypes > 0 ) & 1617 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 1618 & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 1619 1433 1620 END SUBROUTINE dia_obs_dealloc 1434 1621 … … 1436 1623 !! 1437 1624 !! *** ROUTINE ini_date *** 1438 !! 1439 !! ** Purpose : Get initial dat ain double precision YYYYMMDD.HHMMSS format1440 !! 1441 !! ** Method : Get initial dat ain double precision YYYYMMDD.HHMMSS format1442 !! 1443 !! ** Action : Get initial dat ain double precision YYYYMMDD.HHMMSS format1625 !! 1626 !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 1627 !! 1628 !! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format 1629 !! 1630 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1444 1631 !! 1445 1632 !! History : … … 1452 1639 USE phycst, ONLY : & ! Physical constants 1453 1640 & rday 1454 ! USE daymod, ONLY : & ! Time variables1455 ! & nmonth_len1456 1641 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1457 1642 & rdt … … 1460 1645 1461 1646 !! * Arguments 1462 REAL( KIND=dp), INTENT(OUT) :: ddobsini! Initial date in YYYYMMDD.HHMMSS1647 REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 1463 1648 1464 1649 !! * Local declarations … … 1468 1653 INTEGER :: ihou 1469 1654 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 ! !1655 INTEGER :: imday ! Number of days in month. 1656 INTEGER, DIMENSION(12) :: & 1657 & imonth_len ! Length in days of the months of the current year 1658 REAL(wp) :: zdayfrc ! Fraction of day 1659 1660 ! 1661 ! Initial date initialization (year, month, day, hour, minute) 1662 ! (This assumes that the initial date is for 00z)) 1663 ! 1479 1664 iyea = ndate0 / 10000 1480 1665 imon = ( ndate0  iyea * 10000 ) / 100 … … 1483 1668 imin = 0 1484 1669 1485 ! !1486 ! !Compute number of days + number of hours + min since initial time1487 ! !1670 ! 1671 ! Compute number of days + number of hours + min since initial time 1672 ! 1488 1673 iday = iday + ( nit000 1 ) * rdt / rday 1489 1674 zdayfrc = ( nit000 1 ) * rdt / rday … … 1492 1677 imin = int( (zdayfrc * 24  ihou) * 60 ) 1493 1678 1494 ! !1495 ! !Convert number of days (iday) into a real date1496 ! !1679 ! 1680 ! Convert number of days (iday) into a real date 1681 ! 1497 1682 1498 1683 CALL calc_month_len( iyea, imonth_len ) 1499 1684 1500 1685 DO WHILE ( iday > imonth_len(imon) ) 1501 1686 iday = iday  imonth_len(imon) … … 1508 1693 END DO 1509 1694 1510 ! !1511 ! !Convert it into YYYYMMDD.HHMMSS format.1512 ! !1695 ! 1696 ! Convert it into YYYYMMDD.HHMMSS format. 1697 ! 1513 1698 ddobsini = iyea * 10000_dp + imon * 100_dp + & 1514 1699 & iday + ihou * 0.01_dp + imin * 0.0001_dp … … 1520 1705 !! 1521 1706 !! *** ROUTINE fin_date *** 1522 !! 1523 !! ** Purpose : Get final dat ain double precision YYYYMMDD.HHMMSS format1524 !! 1525 !! ** Method : Get final dat ain double precision YYYYMMDD.HHMMSS format1526 !! 1527 !! ** Action : Get final dat ain double precision YYYYMMDD.HHMMSS format1707 !! 1708 !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 1709 !! 1710 !! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format 1711 !! 1712 !! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format 1528 1713 !! 1529 1714 !! History : … … 1535 1720 USE phycst, ONLY : & ! Physical constants 1536 1721 & rday 1537 ! USE daymod, ONLY : & ! Time variables1538 ! & nmonth_len1539 1722 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1540 1723 & rdt … … 1543 1726 1544 1727 !! * Arguments 1545 REAL( KIND=dp), INTENT(OUT) :: ddobsfin! Final date in YYYYMMDD.HHMMSS1728 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1546 1729 1547 1730 !! * Local declarations … … 1551 1734 INTEGER :: ihou 1552 1735 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 1736 INTEGER :: imday ! Number of days in month. 1737 INTEGER, DIMENSION(12) :: & 1738 & imonth_len ! Length in days of the months of the current year 1739 REAL(wp) :: zdayfrc ! Fraction of day 1740 1558 1741 ! 1559 1742 ! Initial date initialization (year, month, day, hour, minute) … … 1565 1748 ihou = 0 1566 1749 imin = 0 1567 1750 1568 1751 ! 1569 1752 ! Compute number of days + number of hours + min since initial time … … 1580 1763 1581 1764 CALL calc_month_len( iyea, imonth_len ) 1582 1765 1583 1766 DO WHILE ( iday > imonth_len(imon) ) 1584 1767 iday = iday  imonth_len(imon) … … 1598 1781 1599 1782 END SUBROUTINE fin_date 1600 1783 1784 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 1785 1786 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1787 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 1788 INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 1789 & ifiles ! Out number of files for each type 1790 CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 1791 & cobstypes ! List of obs types 1792 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 1793 & cfiles ! List of files for all types 1794 1795 !Local variables 1796 INTEGER :: jfile 1797 INTEGER :: jtype 1798 1799 DO jtype = 1, ntypes 1800 1801 ifiles(jtype) = 0 1802 DO jfile = 1, jpmaxnfiles 1803 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 1804 ifiles(jtype) = ifiles(jtype) + 1 1805 END DO 1806 1807 IF ( ifiles(jtype) == 0 ) THEN 1808 CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// & 1809 & ' set to true but no files available to read' ) 1810 ENDIF 1811 1812 IF(lwp) THEN 1813 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1814 DO jfile = 1, ifiles(jtype) 1815 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1816 END DO 1817 ENDIF 1818 1819 END DO 1820 1821 END SUBROUTINE obs_settypefiles 1822 1823 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1824 & n2dint_default, n2dint_type, & 1825 & ravglamscl_type, ravgphiscl_type, & 1826 & lfp_indegs_type, lavnight_type, & 1827 & n2dint, ravglamscl, ravgphiscl, & 1828 & lfpindegs, lavnight ) 1829 1830 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1831 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1832 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1833 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1834 REAL(wp), INTENT(IN) :: & 1835 & ravglamscl_type, & !E/W diameter of obs footprint for this type 1836 & ravgphiscl_type !N/S diameter of obs footprint for this type 1837 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1838 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1839 CHARACTER(len=8), INTENT(IN) :: ctypein 1840 1841 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1842 & n2dint 1843 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 1844 & ravglamscl, ravgphiscl 1845 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 1846 & lfpindegs, lavnight 1847 1848 lavnight(jtype) = lavnight_type 1849 1850 IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 1851 n2dint(jtype) = n2dint_type 1852 ELSE IF ( n2dint_type == 1 ) THEN 1853 n2dint(jtype) = n2dint_default 1854 ELSE 1855 CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 1856 & ' is not available') 1857 ENDIF 1858 1859 ! For averaging observation footprints set options for size of footprint 1860 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 1861 IF ( ravglamscl_type > 0._wp ) THEN 1862 ravglamscl(jtype) = ravglamscl_type 1863 ELSE 1864 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1865 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) 1866 ENDIF 1867 1868 IF ( ravgphiscl_type > 0._wp ) THEN 1869 ravgphiscl(jtype) = ravgphiscl_type 1870 ELSE 1871 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1872 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) 1873 ENDIF 1874 1875 lfpindegs(jtype) = lfp_indegs_type 1876 1877 ENDIF 1878 1879 ! Write out info 1880 IF(lwp) THEN 1881 IF ( n2dint(jtype) <= 4 ) THEN 1882 WRITE(numout,*) ' '//TRIM(ctypein)// & 1883 & ' model counterparts will be interpolated horizontally' 1884 ELSE IF ( n2dint(jtype) <= 6 ) THEN 1885 WRITE(numout,*) ' '//TRIM(ctypein)// & 1886 & ' model counterparts will be averaged horizontally' 1887 WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) 1888 WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) 1889 IF ( lfpindegs(jtype) ) THEN 1890 WRITE(numout,*) ' '//' (in degrees)' 1891 ELSE 1892 WRITE(numout,*) ' '//' (in metres)' 1893 ENDIF 1894 ENDIF 1895 ENDIF 1896 1897 END SUBROUTINE obs_setinterpopts 1898 1601 1899 END MODULE diaobs
Note: See TracChangeset
for help on using the changeset viewer.