Changeset 10120
- Timestamp:
- 2018-09-12T19:12:15+02:00 (4 years ago)
- Location:
- branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM
- Files:
-
- 19 deleted
- 18 edited
- 4 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/CONFIG/SHARED/namelist_ref
r10005 r10120 1224 1224 &namobs ! observation usage switch ('key_diaobs') 1225 1225 !----------------------------------------------------------------------- 1226 ln_t3d = .false. ! Logical switch for T profile observations 1227 ln_s3d = .false. ! Logical switch for S profile observations 1228 ln_ena = .false. ! Logical switch for ENACT insitu data set 1229 ! ! ln_cor Logical switch for Coriolis insitu data set 1230 ln_profb = .false. ! Logical switch for feedback insitu data set 1231 ln_sla = .false. ! Logical switch for SLA observations 1232 1233 ln_sladt = .false. ! Logical switch for AVISO SLA data 1234 1235 ln_slafb = .false. ! Logical switch for feedback SLA data 1236 ! ln_ssh Logical switch for SSH observations 1237 1238 ln_sst = .false. ! Logical switch for SST observations 1239 ln_reysst = .false. ! ln_reysst Logical switch for Reynolds observations 1240 ln_ghrsst = .false. ! ln_ghrsst Logical switch for GHRSST observations 1241 1242 ln_sstfb = .false. ! Logical switch for feedback SST data 1243 ! ln_sss Logical switch for SSS observations 1244 ln_seaice = .false. ! Logical switch for Sea Ice observations 1245 ! ln_vel3d Logical switch for velocity observations 1246 ! ln_velavcur Logical switch for velocity daily av. cur. 1247 ! ln_velhrcur Logical switch for velocity high freq. cur. 1248 ! ln_velavadcp Logical switch for velocity daily av. ADCP 1249 ! ln_velhradcp Logical switch for velocity high freq. ADCP 1250 ! ln_velfb Logical switch for feedback velocity data 1251 ! ln_grid_global Global distribtion of observations 1252 ! ln_grid_search_lookup Logical switch for obs grid search w/lookup table 1253 ! grid_search_file Grid search lookup file header 1254 ! enactfiles ENACT input observation file names 1255 ! coriofiles Coriolis input observation file name 1256 ! ! profbfiles: Profile feedback input observation file name 1257 profbfiles = 'profiles_01.nc' 1258 ! ln_profb_enatim Enact feedback input time setting switch 1259 ! slafilesact Active SLA input observation file name 1260 ! slafilespas Passive SLA input observation file name 1261 ! ! slafbfiles: Feedback SLA input observation file name 1262 slafbfiles = 'sla_01.nc' 1263 ! sstfiles GHRSST input observation file name 1264 ! ! sstfbfiles: Feedback SST input observation file name 1265 sstfbfiles = 'sst_01.nc' 1266 ! seaicefiles Sea Ice input observation file names 1267 seaicefiles = 'seaice_01.nc' 1268 ! velavcurfiles Vel. cur. daily av. input file name 1269 ! velhvcurfiles Vel. cur. high freq. input file name 1270 ! velavadcpfiles Vel. ADCP daily av. input file name 1271 ! velhvadcpfiles Vel. ADCP high freq. input file name 1272 ! velfbfiles Vel. feedback input observation file name 1273 ! dobsini Initial date in window YYYYMMDD.HHMMSS 1274 ! dobsend Final date in window YYYYMMDD.HHMMSS 1275 ! n1dint Type of vertical interpolation method 1276 ! n2dint Type of horizontal interpolation method 1277 ! ln_nea Rejection of observations near land switch 1278 nmsshc = 0 ! MSSH correction scheme 1279 ! mdtcorr MDT correction 1280 ! mdtcutoff MDT cutoff for computed correction 1281 ln_altbias = .false. ! Logical switch for alt bias 1282 ln_ignmis = .true. ! Logical switch for ignoring missing files 1283 ! endailyavtypes ENACT daily average types 1284 ln_grid_global = .true. 1285 ln_grid_search_lookup = .false. 1226 ln_diaobs = .false. ! Logical switch for the observation operator 1227 ln_t3d = .false. ! Logical switch for T profile observations 1228 ln_s3d = .false. ! Logical switch for S profile observations 1229 ln_sla = .false. ! Logical switch for SLA observations 1230 ln_sst = .false. ! Logical switch for SST observations 1231 ln_sic = .false. ! Logical switch for Sea Ice observations 1232 ln_vel3d = .false. ! Logical switch for velocity observations 1233 ln_sss = .false. ! Logical swithc for SSS observations 1234 ln_logchl = .false. ! Logical switch for log(Chl) obs 1235 ln_spm = .false. ! Logical switch for SPM obs 1236 ln_fco2 = .false. 1237 ln_pco2 = .false. 1238 ln_altbias = .false. ! Logical switch for altimeter bias correction 1239 ln_sstbias = .false. ! Logical switch for SST bias correction 1240 ln_nea = .false. ! Logical switch for rejection of observations near land 1241 ln_grid_global = .true. ! Logical switch for global distribution of observations 1242 ln_grid_search_lookup = .false. ! Logical switch for obs grid search w/lookup table 1243 ln_ignmis = .true. ! Logical switch for ignoring missing files 1244 ln_s_at_t = .false. ! Logical switch for computing model S at T obs if not there 1245 ln_sstnight = .false. ! Logical switch for calculating night-time average for SST obs 1246 ln_sla_fp_indegs = .true. ! Logical: T=> averaging footpring is in degrees, F=> in metres 1247 ln_sst_fp_indegs = .true. 1248 ln_sss_fp_indegs = .true. 1249 ln_sic_fp_indegs = .true. 1250 ! All of the *files* variables below are arrays. Use namelist_cfg to add more files 1251 cn_profbfiles = 'profiles_01.nc' ! Profile feedback input observation file names 1252 cn_slafbfiles = 'sla_01.nc' ! SLA feedback input observation file names 1253 cn_sstfbfiles = 'sst_01.nc' ! SST feedback input observation file names 1254 cn_sicfbfiles = 'sic_01.nc' ! SIC feedback input observation file names 1255 cn_velfbfiles = 'vel_01.nc' ! Velocity feedback input observation file names 1256 cn_sssfbfiles = 'sss_01.nc' ! SSS feedback input observation file names 1257 cn_altbiasfile = 'altbias.nc' ! Altimeter bias input file name 1258 cn_sstbiasfiles = 'sstbias.nc' ! SST bias input file names 1259 cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name 1260 rn_gridsearchres = 0.5 ! Grid search resolution 1261 rn_sla_avglamscl = 0. ! E/W diameter of SLA observation footprint (metres/degrees) 1262 rn_sla_avgphiscl = 0. ! N/S diameter of SLA observation footprint (metres/degrees) 1263 rn_sst_avglamscl = 0. 1264 rn_sst_avgphiscl = 0. 1265 rn_sss_avglamscl = 0. 1266 rn_sss_avgphiscl = 0. 1267 rn_sic_avglamscl = 0. 1268 rn_sic_avgphiscl = 0. 1269 nn_1dint = 0 ! Type of vertical interpolation method 1270 nn_2dint = 0 ! Default horizontal interpolation method 1271 nn_2dint_sla = 0 ! Horizontal interpolation method for SLA 1272 nn_2dint_sst = 0 ! Horizontal interpolation method for SST 1273 nn_2dint_sss = 0 ! Horizontal interpolation method for SSS 1274 nn_2dint_sic = 0 ! Horizontal interpolation method for SIC 1275 nn_msshc = 0 ! MSSH correction scheme 1276 rn_mdtcorr = 1.61 ! MDT correction 1277 rn_mdtcutoff = 65.0 ! MDT cutoff for computed correction 1278 nn_profdavtypes = -1 ! Profile daily average types - array 1286 1279 / 1287 1280 !----------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90
r8447 r10120 243 243 ! Set additional default values (note that most values are set in the reference namelist) 244 244 245 ! This doesn't work because nitiaufin hasn't been set yet 245 246 IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin 246 247 … … 318 319 IF ( (.NOT. ln_incpc) .AND. ln_incpc_only) & 319 320 & CALL ctl_stop (' if you set ln_incpc_only to .true. then you need to set ln_incpc to .true. as well' ) 320 321 WRITE(numout,*) ' time step is = ',rdt,'you choose to write pcbias at nn_bias_itwrt = ',nn_bias_itwrt,'and end of iau is rday/rdt=',rday/rdt 321 322 ! This doesn't work because nitiaufin hasn't been set yet (but the old code with rday/rdt was also wrong) 323 ! WRITE(numout,*) ' time step is = ',rdt,'you choose to write pcbias at nn_bias_itwrt = ',nn_bias_itwrt,'and end of iau is nitiaufin=',nitiaufin 322 324 ENDIF 323 325 IF( .NOT. ln_bias ) RETURN -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r10005 r10120 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 … … 53 44 & calc_date ! Compute the date of a timestep 54 45 55 !! * Shared Module variables56 LOGICAL, PUBLIC, PARAMETER :: &57 #if defined key_diaobs58 & lk_diaobs = .TRUE. !: Logical switch for observation diangostics59 #else60 & lk_diaobs = .FALSE. !: Logical switch for observation diangostics61 #endif62 63 46 !! * Module variables 64 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles 65 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles 66 LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set 67 LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set 68 LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles 69 LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies 70 LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files 71 LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files 72 LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature 73 LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature 74 LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data 75 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files 76 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration 77 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations 78 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data 79 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data 80 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data 81 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data 82 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files 83 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height 84 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity 85 LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations 86 LOGICAL, PUBLIC :: ln_nea !: Remove observations near land 87 LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias 88 LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files 89 LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations 90 91 REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS 92 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS 93 94 INTEGER, PUBLIC :: n1dint !: Vertical interpolation method 95 INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method 96 47 LOGICAL, PUBLIC :: & 48 & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. 49 LOGICAL :: ln_diaobs !: Logical switch for the obs operator 50 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 51 LOGICAL :: ln_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=> sea-ice obs footprint size specified in degrees, F=> in metres 55 56 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 57 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 58 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 59 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 60 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 61 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 62 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 63 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 64 65 INTEGER :: nn_1dint !: Vertical interpolation method 66 INTEGER :: nn_2dint !: Default horizontal interpolation method 67 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method 68 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method 69 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method 70 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method 71 97 72 INTEGER, DIMENSION(imaxavtypes) :: & 98 & endailyavtypes !: ENACT data types which are daily average 99 100 INTEGER, PARAMETER :: MaxNumFiles = 1000 101 LOGICAL, DIMENSION(MaxNumFiles) :: & 102 & ln_profb_ena, & !: Is the feedback files from ENACT data ? 103 ! !: If so use endailyavtypes 104 & ln_profb_enatim !: Change tim for 820 enact data set. 105 106 LOGICAL, DIMENSION(MaxNumFiles) :: & 107 & ln_velfb_av !: Is the velocity feedback files daily average? 73 & nn_profdavtypes !: Profile data types representing a daily average 74 INTEGER :: nproftypes !: Number of profile obs types 75 INTEGER :: nsurftypes !: Number of surface obs types 76 INTEGER, DIMENSION(:), ALLOCATABLE :: & 77 & nvarsprof, & !: Number of profile variables 78 & nvarssurf !: Number of surface variables 79 INTEGER, DIMENSION(:), ALLOCATABLE :: & 80 & nextrprof, & !: Number of profile extra variables 81 & nextrsurf !: Number of surface extra variables 82 INTEGER, DIMENSION(:), ALLOCATABLE :: & 83 & n2dintsurf !: Interpolation option for surface variables 84 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 85 & ravglamscl, & !: E/W diameter of averaging footprint for surface variables 86 & ravgphiscl !: N/S diameter of averaging footprint for surface variables 108 87 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 109 & ld_enact !: Profile data is ENACT so use endailyavtypes 110 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 111 & ld_velav !: Velocity data is daily averaged 112 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 113 & ld_sstnight !: SST observation corresponds to night mean 88 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 89 & llnightav !: Logical for calculating night-time averages 90 91 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 92 & surfdata, & !: Initial surface data 93 & surfdataqc !: Surface data after quality control 94 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 95 & profdata, & !: Initial profile data 96 & profdataqc !: Profile data after quality control 97 98 CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 99 & cobstypesprof, & !: Profile obs types 100 & cobstypessurf !: Surface obs types 114 101 115 102 !!---------------------------------------------------------------------- … … 119 106 !!---------------------------------------------------------------------- 120 107 108 !! * Substitutions 109 # include "domzgr_substitute.h90" 121 110 CONTAINS 122 111 … … 136 125 !! ! 06-10 (A. Weaver) Cleaning and add controls 137 126 !! ! 07-03 (K. Mogensen) General handling of profiles 127 !! ! 14-08 (J.While) Incorporated SST bias correction 128 !! ! 15-02 (M. Martin) Simplification of namelist and code 138 129 !!---------------------------------------------------------------------- 139 130 … … 141 132 142 133 !! * Local declarations 143 CHARACTER(len=128) :: enactfiles(MaxNumFiles) 144 CHARACTER(len=128) :: coriofiles(MaxNumFiles) 145 CHARACTER(len=128) :: profbfiles(MaxNumFiles) 146 CHARACTER(len=128) :: sstfiles(MaxNumFiles) 147 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 148 CHARACTER(len=128) :: slafilesact(MaxNumFiles) 149 CHARACTER(len=128) :: slafilespas(MaxNumFiles) 150 CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 151 CHARACTER(len=128) :: seaicefiles(MaxNumFiles) 152 CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 153 CHARACTER(len=128) :: veladcpfiles(MaxNumFiles) 154 CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 155 CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 156 CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 157 CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 158 CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 159 CHARACTER(LEN=128) :: reysstname 160 CHARACTER(LEN=12) :: reysstfmt 161 CHARACTER(LEN=128) :: bias_file 162 CHARACTER(LEN=20) :: datestr=" ", timestr=" " 163 NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, & 164 & ln_sla, ln_sladt, ln_slafb, & 165 & ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, & 166 & enactfiles, coriofiles, profbfiles, & 167 & slafilesact, slafilespas, slafbfiles, & 168 & sstfiles, sstfbfiles, & 169 & ln_seaice, seaicefiles, & 170 & dobsini, dobsend, n1dint, n2dint, & 171 & nmsshc, mdtcorr, mdtcutoff, & 172 & ln_reysst, ln_ghrsst, reysstname, reysstfmt, & 134 INTEGER, PARAMETER :: & 135 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 136 INTEGER, DIMENSION(:), ALLOCATABLE :: & 137 & ifilesprof, & ! Number of profile files 138 & ifilessurf ! Number of surface files 139 INTEGER :: ios ! Local integer output status for namelist read 140 INTEGER :: jtype ! Counter for obs types 141 INTEGER :: jvar ! Counter for variables 142 INTEGER :: jfile ! Counter for files 143 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 144 145 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 146 & cn_profbfiles, & ! T/S profile input filenames 147 & cn_sstfbfiles, & ! Sea surface temperature input filenames 148 & cn_slafbfiles, & ! Sea level anomaly input filenames 149 & cn_sicfbfiles, & ! Seaice concentration input filenames 150 & cn_velfbfiles, & ! Velocity profile input filenames 151 & cn_sssfbfiles, & ! Sea surface salinity input filenames 152 & cn_logchlfbfiles, & ! Log(Chl) input filenames 153 & cn_spmfbfiles, & ! Sediment input filenames 154 & cn_fco2fbfiles, & ! fco2 input filenames 155 & cn_pco2fbfiles, & ! pco2 input filenames 156 & cn_sstbiasfiles ! SST bias input filenames 157 158 CHARACTER(LEN=128) :: & 159 & cn_altbiasfile ! Altimeter bias input filename 160 161 162 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 163 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 164 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 165 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 166 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 167 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 168 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 169 LOGICAL :: ln_logchl ! Logical switch for log(Chl) obs 170 LOGICAL :: ln_spm ! Logical switch for sediment obs 171 LOGICAL :: ln_fco2 ! Logical switch for fco2 obs 172 LOGICAL :: ln_pco2 ! Logical switch for pco2 obs 173 LOGICAL :: ln_nea ! Logical switch to remove obs near land 174 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 175 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 176 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 177 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 178 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 179 180 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 181 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 182 183 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 184 & clproffiles, & ! Profile filenames 185 & clsurffiles ! Surface filenames 186 187 LOGICAL :: llvar1 ! Logical for profile variable 1 188 LOGICAL :: llvar2 ! Logical for profile variable 1 189 190 REAL(wp), POINTER, DIMENSION(:,:) :: & 191 & zglam1, & ! Model longitudes for profile variable 1 192 & zglam2 ! Model longitudes for profile variable 2 193 REAL(wp), POINTER, DIMENSION(:,:) :: & 194 & zgphi1, & ! Model latitudes for profile variable 1 195 & zgphi2 ! Model latitudes for profile variable 2 196 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 197 & zmask1, & ! Model land/sea mask associated with variable 1 198 & zmask2 ! Model land/sea mask associated with variable 2 199 200 201 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 202 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 203 & ln_logchl, ln_spm, ln_fco2, ln_pco2, & 204 & ln_altbias, ln_sstbias, ln_nea, & 205 & ln_grid_global, ln_grid_search_lookup, & 206 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 173 207 & ln_sstnight, & 174 & ln_grid_search_lookup, & 175 & grid_search_file, grid_search_res, & 176 & ln_grid_global, bias_file, ln_altbias, & 177 & endailyavtypes, ln_s_at_t, ln_profb_ena, & 178 & ln_vel3d, ln_velavcur, velavcurfiles, & 179 & ln_velhrcur, velhrcurfiles, & 180 & ln_velavadcp, velavadcpfiles, & 181 & ln_velhradcp, velhradcpfiles, & 182 & ln_velfb, velfbfiles, ln_velfb_av, & 183 & ln_profb_enatim, ln_ignmis, ln_cl4 184 185 INTEGER :: jprofset 186 INTEGER :: jveloset 187 INTEGER :: jvar 188 INTEGER :: jnumenact 189 INTEGER :: jnumcorio 190 INTEGER :: jnumprofb 191 INTEGER :: jnumslaact 192 INTEGER :: jnumslapas 193 INTEGER :: jnumslafb 194 INTEGER :: jnumsst 195 INTEGER :: jnumsstfb 196 INTEGER :: jnumseaice 197 INTEGER :: jnumvelavcur 198 INTEGER :: jnumvelhrcur 199 INTEGER :: jnumvelavadcp 200 INTEGER :: jnumvelhradcp 201 INTEGER :: jnumvelfb 202 INTEGER :: ji 203 INTEGER :: jset 204 INTEGER :: ios ! Local integer output status for namelist read 205 LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 208 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 209 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 210 & cn_profbfiles, cn_slafbfiles, & 211 & cn_sstfbfiles, cn_sicfbfiles, & 212 & cn_velfbfiles, cn_sssfbfiles, & 213 & cn_logchlfbfiles, cn_spmfbfiles, & 214 & cn_fco2fbfiles, cn_pco2fbfiles, & 215 & cn_sstbiasfiles, cn_altbiasfile, & 216 & cn_gridsearchfile, rn_gridsearchres, & 217 & rn_dobsini, rn_dobsend, & 218 & rn_sla_avglamscl, rn_sla_avgphiscl, & 219 & rn_sst_avglamscl, rn_sst_avgphiscl, & 220 & rn_sss_avglamscl, rn_sss_avgphiscl, & 221 & rn_sic_avglamscl, rn_sic_avgphiscl, & 222 & nn_1dint, nn_2dint, & 223 & nn_2dint_sla, nn_2dint_sst, & 224 & nn_2dint_sss, nn_2dint_sic, & 225 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 226 & nn_profdavtypes 227 228 CALL wrk_alloc( jpi, jpj, zglam1 ) 229 CALL wrk_alloc( jpi, jpj, zglam2 ) 230 CALL wrk_alloc( jpi, jpj, zgphi1 ) 231 CALL wrk_alloc( jpi, jpj, zgphi2 ) 232 CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 233 CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 206 234 207 235 !----------------------------------------------------------------------- … … 209 237 !----------------------------------------------------------------------- 210 238 211 enactfiles(:) = '' 212 coriofiles(:) = '' 213 profbfiles(:) = '' 214 slafilesact(:) = '' 215 slafilespas(:) = '' 216 slafbfiles(:) = '' 217 sstfiles(:) = '' 218 sstfbfiles(:) = '' 219 seaicefiles(:) = '' 220 velcurfiles(:) = '' 221 veladcpfiles(:) = '' 222 velavcurfiles(:) = '' 223 velhrcurfiles(:) = '' 224 velavadcpfiles(:) = '' 225 velhradcpfiles(:) = '' 226 velfbfiles(:) = '' 227 velcurfiles(:) = '' 228 veladcpfiles(:) = '' 229 endailyavtypes(:) = -1 230 endailyavtypes(1) = 820 231 ln_profb_ena(:) = .FALSE. 232 ln_profb_enatim(:) = .TRUE. 233 ln_velfb_av(:) = .FALSE. 234 ln_ignmis = .FALSE. 235 236 CALL ini_date( dobsini ) 237 CALL fin_date( dobsend ) 238 239 ! Read Namelist namobs : control observation diagnostics 240 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation 239 ! Some namelist arrays need initialising 240 cn_profbfiles(:) = '' 241 cn_slafbfiles(:) = '' 242 cn_sstfbfiles(:) = '' 243 cn_sicfbfiles(:) = '' 244 cn_velfbfiles(:) = '' 245 cn_sssfbfiles(:) = '' 246 cn_logchlfbfiles(:) = '' 247 cn_spmfbfiles(:) = '' 248 cn_fco2fbfiles(:) = '' 249 cn_pco2fbfiles(:) = '' 250 cn_sstbiasfiles(:) = '' 251 nn_profdavtypes(:) = -1 252 253 CALL ini_date( rn_dobsini ) 254 CALL fin_date( rn_dobsend ) 255 256 ! Read namelist namobs : control observation diagnostics 257 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 241 258 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 242 259 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 243 260 244 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation261 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 245 262 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 246 263 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 247 264 IF(lwm) WRITE ( numond, namobs ) 248 265 249 ! Count number of files for each type 250 IF (ln_ena) THEN 251 lmask(:) = .FALSE. 252 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 253 jnumenact = COUNT(lmask) 254 ENDIF 255 IF (ln_cor) THEN 256 lmask(:) = .FALSE. 257 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 258 jnumcorio = COUNT(lmask) 259 ENDIF 260 IF (ln_profb) THEN 261 lmask(:) = .FALSE. 262 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 263 jnumprofb = COUNT(lmask) 264 ENDIF 265 IF (ln_sladt) THEN 266 lmask(:) = .FALSE. 267 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 268 jnumslaact = COUNT(lmask) 269 lmask(:) = .FALSE. 270 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 271 jnumslapas = COUNT(lmask) 272 ENDIF 273 IF (ln_slafb) THEN 274 lmask(:) = .FALSE. 275 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 276 jnumslafb = COUNT(lmask) 277 lmask(:) = .FALSE. 278 ENDIF 279 IF (ln_ghrsst) THEN 280 lmask(:) = .FALSE. 281 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 282 jnumsst = COUNT(lmask) 283 ENDIF 284 IF (ln_sstfb) THEN 285 lmask(:) = .FALSE. 286 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 287 jnumsstfb = COUNT(lmask) 288 lmask(:) = .FALSE. 289 ENDIF 290 IF (ln_seaice) THEN 291 lmask(:) = .FALSE. 292 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 293 jnumseaice = COUNT(lmask) 294 ENDIF 295 IF (ln_velavcur) THEN 296 lmask(:) = .FALSE. 297 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 298 jnumvelavcur = COUNT(lmask) 299 ENDIF 300 IF (ln_velhrcur) THEN 301 lmask(:) = .FALSE. 302 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 303 jnumvelhrcur = COUNT(lmask) 304 ENDIF 305 IF (ln_velavadcp) THEN 306 lmask(:) = .FALSE. 307 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 308 jnumvelavadcp = COUNT(lmask) 309 ENDIF 310 IF (ln_velhradcp) THEN 311 lmask(:) = .FALSE. 312 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 313 jnumvelhradcp = COUNT(lmask) 314 ENDIF 315 IF (ln_velfb) THEN 316 lmask(:) = .FALSE. 317 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 318 jnumvelfb = COUNT(lmask) 319 lmask(:) = .FALSE. 320 ENDIF 321 322 ! Control print 266 lk_diaobs = .FALSE. 267 #if defined key_diaobs 268 IF ( ln_diaobs ) lk_diaobs = .TRUE. 269 #endif 270 271 IF ( .NOT. lk_diaobs ) THEN 272 IF(lwp) WRITE(numout,cform_war) 273 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' 274 RETURN 275 ENDIF 276 323 277 IF(lwp) THEN 324 278 WRITE(numout,*) … … 326 280 WRITE(numout,*) '~~~~~~~~~~~~' 327 281 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 328 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 329 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 330 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena 331 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor 332 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb 333 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 334 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt 335 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb 336 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 337 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 338 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst 339 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst 340 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb 341 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 342 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 343 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 344 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 345 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 346 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 347 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 348 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 349 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 350 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 351 WRITE(numout,*) & 352 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 282 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 283 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 284 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 285 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 286 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 287 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 288 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 289 WRITE(numout,*) ' Logical switch for log(Chl) observations ln_logchl = ', ln_logchl 290 WRITE(numout,*) ' Logical switch for SPM observations ln_spm = ', ln_spm 291 WRITE(numout,*) ' Logical switch for FCO2 observations ln_fco2 = ', ln_fco2 292 WRITE(numout,*) ' Logical switch for PCO2 observations ln_pco2 = ', ln_pco2 293 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 294 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 353 295 IF (ln_grid_search_lookup) & 354 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 355 IF (ln_ena) THEN 356 DO ji = 1, jnumenact 357 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 358 TRIM(enactfiles(ji)) 359 END DO 296 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 297 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 298 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 299 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 300 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 301 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 302 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 303 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 304 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 305 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 306 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 307 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 308 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 309 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 310 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 311 ENDIF 312 !----------------------------------------------------------------------- 313 ! Set up list of observation types to be used 314 ! and the files associated with each type 315 !----------------------------------------------------------------------- 316 317 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 318 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 319 & ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) 320 321 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 322 IF(lwp) WRITE(numout,cform_war) 323 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 324 & ' are set to .FALSE. so turning off calls to dia_obs' 325 nwarn = nwarn + 1 326 lk_diaobs = .FALSE. 327 RETURN 328 ENDIF 329 330 IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes 331 IF ( nproftypes > 0 ) THEN 332 333 ALLOCATE( cobstypesprof(nproftypes) ) 334 ALLOCATE( ifilesprof(nproftypes) ) 335 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 336 337 jtype = 0 338 IF (ln_t3d .OR. ln_s3d) THEN 339 jtype = jtype + 1 340 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', & 341 & cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 360 342 ENDIF 361 IF (ln_cor) THEN 362 DO ji = 1, jnumcorio 363 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 364 TRIM(coriofiles(ji)) 365 END DO 343 IF (ln_vel3d) THEN 344 jtype = jtype + 1 345 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', & 346 & cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 366 347 ENDIF 367 IF (ln_profb) THEN 368 DO ji = 1, jnumprofb 369 IF (ln_profb_ena(ji)) THEN 370 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 371 TRIM(profbfiles(ji)) 372 ELSE 373 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 374 TRIM(profbfiles(ji)) 375 ENDIF 376 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 377 END DO 348 349 ENDIF 350 351 IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes 352 IF ( nsurftypes > 0 ) THEN 353 354 ALLOCATE( cobstypessurf(nsurftypes) ) 355 ALLOCATE( ifilessurf(nsurftypes) ) 356 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 357 ALLOCATE(n2dintsurf(nsurftypes)) 358 ALLOCATE(ravglamscl(nsurftypes)) 359 ALLOCATE(ravgphiscl(nsurftypes)) 360 ALLOCATE(lfpindegs(nsurftypes)) 361 ALLOCATE(llnightav(nsurftypes)) 362 363 jtype = 0 364 IF (ln_sla) THEN 365 jtype = jtype + 1 366 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', & 367 & cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 368 CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', & 369 & nn_2dint, nn_2dint_sla, & 370 & rn_sla_avglamscl, rn_sla_avgphiscl, & 371 & ln_sla_fp_indegs, .FALSE., & 372 & n2dintsurf, ravglamscl, ravgphiscl, & 373 & lfpindegs, llnightav ) 378 374 ENDIF 379 IF (ln_sladt) THEN 380 DO ji = 1, jnumslaact 381 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 382 TRIM(slafilesact(ji)) 383 END DO 384 DO ji = 1, jnumslapas 385 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 386 TRIM(slafilespas(ji)) 387 END DO 375 IF (ln_sst) THEN 376 jtype = jtype + 1 377 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', & 378 & cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 379 CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', & 380 & nn_2dint, nn_2dint_sst, & 381 & rn_sst_avglamscl, rn_sst_avgphiscl, & 382 & ln_sst_fp_indegs, ln_sstnight, & 383 & n2dintsurf, ravglamscl, ravgphiscl, & 384 & lfpindegs, llnightav ) 388 385 ENDIF 389 IF (ln_slafb) THEN 390 DO ji = 1, jnumslafb 391 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 392 TRIM(slafbfiles(ji)) 393 END DO 386 #if defined key_lim2 || defined key_lim3 || defined key_cice 387 IF (ln_sic) THEN 388 jtype = jtype + 1 389 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', & 390 & cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 391 CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', & 392 & nn_2dint, nn_2dint_sic, & 393 & rn_sic_avglamscl, rn_sic_avgphiscl, & 394 & ln_sic_fp_indegs, .FALSE., & 395 & n2dintsurf, ravglamscl, ravgphiscl, & 396 & lfpindegs, llnightav ) 394 397 ENDIF 395 IF (ln_ghrsst) THEN 396 DO ji = 1, jnumsst 397 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 398 TRIM(sstfiles(ji)) 399 END DO 398 #endif 399 IF (ln_sss) THEN 400 jtype = jtype + 1 401 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', & 402 & cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 403 CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', & 404 & nn_2dint, nn_2dint_sss, & 405 & rn_sss_avglamscl, rn_sss_avgphiscl, & 406 & ln_sss_fp_indegs, .FALSE., & 407 & n2dintsurf, ravglamscl, ravgphiscl, & 408 & lfpindegs, llnightav ) 400 409 ENDIF 401 IF (ln_sstfb) THEN 402 DO ji = 1, jnumsstfb 403 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 404 TRIM(sstfbfiles(ji)) 405 END DO 410 411 IF (ln_logchl) THEN 412 jtype = jtype + 1 413 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', & 414 & cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 415 CALL obs_setinterpopts( nsurftypes, jtype, 'logchl', & 416 & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 417 & n2dintsurf, ravglamscl, ravgphiscl, & 418 & lfpindegs, llnightav ) 406 419 ENDIF 407 IF (ln_seaice) THEN 408 DO ji = 1, jnumseaice 409 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 410 TRIM(seaicefiles(ji)) 411 END DO 420 421 IF (ln_spm) THEN 422 jtype = jtype + 1 423 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm ', & 424 & cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 425 CALL obs_setinterpopts( nsurftypes, jtype, 'spm ', & 426 & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 427 & n2dintsurf, ravglamscl, ravgphiscl, & 428 & lfpindegs, llnightav ) 412 429 ENDIF 413 IF (ln_velavcur) THEN 414 DO ji = 1, jnumvelavcur 415 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 416 TRIM(velavcurfiles(ji)) 417 END DO 430 431 IF (ln_fco2) THEN 432 jtype = jtype + 1 433 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2 ', & 434 & cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 435 CALL obs_setinterpopts( nsurftypes, jtype, 'fco2 ', & 436 & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 437 & n2dintsurf, ravglamscl, ravgphiscl, & 438 & lfpindegs, llnightav ) 418 439 ENDIF 419 IF (ln_velhrcur) THEN 420 DO ji = 1, jnumvelhrcur 421 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 422 TRIM(velhrcurfiles(ji)) 423 END DO 440 441 IF (ln_pco2) THEN 442 jtype = jtype + 1 443 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2 ', & 444 & cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 445 CALL obs_setinterpopts( nsurftypes, jtype, 'pco2 ', & 446 & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 447 & n2dintsurf, ravglamscl, ravgphiscl, & 448 & lfpindegs, llnightav ) 424 449 ENDIF 425 IF (ln_velavadcp) THEN 426 DO ji = 1, jnumvelavadcp 427 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 428 TRIM(velavadcpfiles(ji)) 429 END DO 430 ENDIF 431 IF (ln_velhradcp) THEN 432 DO ji = 1, jnumvelhradcp 433 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 434 TRIM(velhradcpfiles(ji)) 435 END DO 436 ENDIF 437 IF (ln_velfb) THEN 438 DO ji = 1, jnumvelfb 439 IF (ln_velfb_av(ji)) THEN 440 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 441 TRIM(velfbfiles(ji)) 442 ELSE 443 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 444 TRIM(velfbfiles(ji)) 445 ENDIF 446 END DO 447 ENDIF 448 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 449 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 450 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint 451 WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint 452 WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea 453 WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc 454 WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr 455 WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff 456 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 457 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 458 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 459 460 ENDIF 461 450 451 ENDIF 452 453 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 454 455 456 !----------------------------------------------------------------------- 457 ! Obs operator parameter checking and initialisations 458 !----------------------------------------------------------------------- 459 462 460 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 463 461 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 465 463 ENDIF 466 464 467 CALL obs_typ_init 468 469 CALL mppmap_init 470 471 ! Parameter control 472 #if defined key_diaobs 473 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 474 & ( .NOT. ln_vel3d ).AND. & 475 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 476 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 477 IF(lwp) WRITE(numout,cform_war) 478 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 479 & ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 480 nwarn = nwarn + 1 481 ENDIF 482 #endif 483 484 CALL obs_grid_setup( ) 485 IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 465 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 486 466 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 487 467 & ' is not available') 488 468 ENDIF 489 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 469 470 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 490 471 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 491 472 & ' is not available') 492 473 ENDIF 474 475 CALL obs_typ_init 476 477 CALL mppmap_init 478 479 CALL obs_grid_setup( ) 493 480 494 481 !----------------------------------------------------------------------- 495 482 ! Depending on switches read the various observation types 496 483 !----------------------------------------------------------------------- 497 ! - Temperature/salinity profiles 498 499 IF ( ln_t3d .OR. ln_s3d ) THEN 500 501 ! Set the number of variables for profiles to 2 (T and S) 502 nprofvars = 2 503 ! Set the number of extra variables for profiles to 1 (insitu temp). 504 nprofextr = 1 505 506 ! Count how may insitu data sets we have and allocate data. 507 jprofset = 0 508 IF ( ln_ena ) jprofset = jprofset + 1 509 IF ( ln_cor ) jprofset = jprofset + 1 510 IF ( ln_profb ) jprofset = jprofset + jnumprofb 511 nprofsets = jprofset 512 IF ( nprofsets > 0 ) THEN 513 ALLOCATE(ld_enact(nprofsets)) 514 ALLOCATE(profdata(nprofsets)) 515 ALLOCATE(prodatqc(nprofsets)) 516 ENDIF 517 518 jprofset = 0 519 520 ! ENACT insitu data 521 522 IF ( ln_ena ) THEN 523 524 jprofset = jprofset + 1 525 526 ld_enact(jprofset) = .TRUE. 527 528 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 529 & jnumenact, enactfiles(1:jnumenact), & 530 & nprofvars, nprofextr, & 531 & nitend-nit000+2, & 532 & dobsini, dobsend, ln_t3d, ln_s3d, & 533 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 534 & kdailyavtypes = endailyavtypes ) 535 536 DO jvar = 1, 2 537 538 CALL obs_prof_staend( profdata(jprofset), jvar ) 539 484 485 IF ( nproftypes > 0 ) THEN 486 487 ALLOCATE(profdata(nproftypes)) 488 ALLOCATE(profdataqc(nproftypes)) 489 ALLOCATE(nvarsprof(nproftypes)) 490 ALLOCATE(nextrprof(nproftypes)) 491 492 DO jtype = 1, nproftypes 493 494 nvarsprof(jtype) = 2 495 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 496 nextrprof(jtype) = 1 497 llvar1 = ln_t3d 498 llvar2 = ln_s3d 499 zglam1 = glamt 500 zgphi1 = gphit 501 zmask1 = tmask 502 zglam2 = glamt 503 zgphi2 = gphit 504 zmask2 = tmask 505 ENDIF 506 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 507 nextrprof(jtype) = 2 508 llvar1 = ln_vel3d 509 llvar2 = ln_vel3d 510 zglam1 = glamu 511 zgphi1 = gphiu 512 zmask1 = umask 513 zglam2 = glamv 514 zgphi2 = gphiv 515 zmask2 = vmask 516 ENDIF 517 518 !Read in profile or profile obs types 519 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 520 & clproffiles(jtype,1:ifilesprof(jtype)), & 521 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 522 & rn_dobsini, rn_dobsend, llvar1, llvar2, & 523 & ln_ignmis, ln_s_at_t, .FALSE., & 524 & kdailyavtypes = nn_profdavtypes ) 525 526 DO jvar = 1, nvarsprof(jtype) 527 CALL obs_prof_staend( profdata(jtype), jvar ) 540 528 END DO 541 529 542 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 543 & ln_t3d, ln_s3d, ln_nea, & 544 & kdailyavtypes=endailyavtypes ) 545 546 ENDIF 547 548 ! Coriolis insitu data 549 550 IF ( ln_cor ) THEN 551 552 jprofset = jprofset + 1 553 554 ld_enact(jprofset) = .FALSE. 555 556 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 557 & jnumcorio, coriofiles(1:jnumcorio), & 558 & nprofvars, nprofextr, & 559 & nitend-nit000+2, & 560 & dobsini, dobsend, ln_t3d, ln_s3d, & 561 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 562 563 DO jvar = 1, 2 564 565 CALL obs_prof_staend( profdata(jprofset), jvar ) 566 567 END DO 568 569 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 570 & ln_t3d, ln_s3d, ln_nea ) 571 572 ENDIF 573 574 ! Feedback insitu data 575 576 IF ( ln_profb ) THEN 577 578 DO jset = 1, jnumprofb 579 580 jprofset = jprofset + 1 581 ld_enact (jprofset) = ln_profb_ena(jset) 582 583 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 584 & 1, profbfiles(jset:jset), & 585 & nprofvars, nprofextr, & 586 & nitend-nit000+2, & 587 & dobsini, dobsend, ln_t3d, ln_s3d, & 588 & ln_ignmis, ln_s_at_t, & 589 & ld_enact(jprofset).AND.& 590 & ln_profb_enatim(jset), & 591 & .FALSE., kdailyavtypes = endailyavtypes ) 592 593 DO jvar = 1, 2 594 595 CALL obs_prof_staend( profdata(jprofset), jvar ) 596 530 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 531 & llvar1, llvar2, & 532 & jpi, jpj, jpk, & 533 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 534 & ln_nea, ln_bound_reject, & 535 & kdailyavtypes = nn_profdavtypes ) 536 537 END DO 538 539 DEALLOCATE( ifilesprof, clproffiles ) 540 541 ENDIF 542 543 IF ( nsurftypes > 0 ) THEN 544 545 ALLOCATE(surfdata(nsurftypes)) 546 ALLOCATE(surfdataqc(nsurftypes)) 547 ALLOCATE(nvarssurf(nsurftypes)) 548 ALLOCATE(nextrsurf(nsurftypes)) 549 550 DO jtype = 1, nsurftypes 551 552 nvarssurf(jtype) = 1 553 nextrsurf(jtype) = 0 554 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 555 556 !Read in surface obs types 557 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 558 & clsurffiles(jtype,1:ifilessurf(jtype)), & 559 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 560 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 561 562 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 563 564 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 565 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 566 IF ( ln_altbias ) & 567 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 568 ENDIF 569 570 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 571 jnumsstbias = 0 572 DO jfile = 1, jpmaxnfiles 573 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 574 & jnumsstbias = jnumsstbias + 1 597 575 END DO 598 599 IF ( ld_enact(jprofset) ) THEN 600 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 601 & ln_t3d, ln_s3d, ln_nea, & 602 & kdailyavtypes = endailyavtypes ) 603 ELSE 604 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 605 & ln_t3d, ln_s3d, ln_nea ) 576 IF ( jnumsstbias == 0 ) THEN 577 CALL ctl_stop("ln_sstbias set but no bias files to read in") 606 578 ENDIF 607 608 END DO 609 610 ENDIF 611 612 ENDIF 613 614 ! - Sea level anomalies 615 IF ( ln_sla ) THEN 616 ! Set the number of variables for sla to 1 617 nslavars = 1 618 619 ! Set the number of extra variables for sla to 2 620 nslaextr = 2 621 622 ! Set the number of sla data sets to 2 623 nslasets = 0 624 IF ( ln_sladt ) THEN 625 nslasets = nslasets + 2 626 ENDIF 627 IF ( ln_slafb ) THEN 628 nslasets = nslasets + jnumslafb 629 ENDIF 630 631 ALLOCATE(sladata(nslasets)) 632 ALLOCATE(sladatqc(nslasets)) 633 sladata(:)%nsurf=0 634 sladatqc(:)%nsurf=0 635 636 nslasets = 0 637 638 ! AVISO SLA data 639 640 IF ( ln_sladt ) THEN 641 642 ! Active SLA observations 643 644 nslasets = nslasets + 1 645 646 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 647 & slafilesact(1:jnumslaact), & 648 & nslavars, nslaextr, nitend-nit000+2, & 649 & dobsini, dobsend, ln_ignmis, .FALSE. ) 650 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 651 & ln_sla, ln_nea ) 652 653 ! Passive SLA observations 654 655 nslasets = nslasets + 1 656 657 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 658 & slafilespas(1:jnumslapas), & 659 & nslavars, nslaextr, nitend-nit000+2, & 660 & dobsini, dobsend, ln_ignmis, .FALSE. ) 661 662 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 663 & ln_sla, ln_nea ) 664 665 ENDIF 666 667 ! Feedback SLA data 668 669 IF ( ln_slafb ) THEN 670 671 DO jset = 1, jnumslafb 672 673 nslasets = nslasets + 1 674 675 CALL obs_rea_sla( 0, sladata(nslasets), 1, & 676 & slafbfiles(jset:jset), & 677 & nslavars, nslaextr, nitend-nit000+2, & 678 & dobsini, dobsend, ln_ignmis, .FALSE. ) 679 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 680 & ln_sla, ln_nea ) 681 682 END DO 683 684 ENDIF 685 686 CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 687 688 ! read in altimeter bias 689 690 IF ( ln_altbias ) THEN 691 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 692 ENDIF 693 694 ENDIF 695 696 ! - Sea surface height 697 IF ( ln_ssh ) THEN 698 IF(lwp) WRITE(numout,*) ' SSH currently not available' 699 ENDIF 700 701 ! - Sea surface temperature 702 IF ( ln_sst ) THEN 703 704 ! Set the number of variables for sst to 1 705 nsstvars = 1 706 707 ! Set the number of extra variables for sst to 0 708 nsstextr = 0 709 710 nsstsets = 0 711 712 IF (ln_reysst) nsstsets = nsstsets + 1 713 IF (ln_ghrsst) nsstsets = nsstsets + 1 714 IF ( ln_sstfb ) THEN 715 nsstsets = nsstsets + jnumsstfb 716 ENDIF 717 718 ALLOCATE(sstdata(nsstsets)) 719 ALLOCATE(sstdatqc(nsstsets)) 720 ALLOCATE(ld_sstnight(nsstsets)) 721 sstdata(:)%nsurf=0 722 sstdatqc(:)%nsurf=0 723 ld_sstnight(:)=.false. 724 725 nsstsets = 0 726 727 IF (ln_reysst) THEN 728 729 nsstsets = nsstsets + 1 730 731 ld_sstnight(nsstsets) = ln_sstnight 732 733 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 734 & nsstvars, nsstextr, & 735 & nitend-nit000+2, dobsini, dobsend ) 736 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 737 & ln_nea ) 738 739 ENDIF 740 741 IF (ln_ghrsst) THEN 742 743 nsstsets = nsstsets + 1 744 745 ld_sstnight(nsstsets) = ln_sstnight 746 747 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 748 & sstfiles(1:jnumsst), & 749 & nsstvars, nsstextr, nitend-nit000+2, & 750 & dobsini, dobsend, ln_ignmis, .FALSE. ) 751 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 752 & ln_nea ) 753 754 ENDIF 755 756 ! Feedback SST data 757 758 IF ( ln_sstfb ) THEN 759 760 DO jset = 1, jnumsstfb 761 762 nsstsets = nsstsets + 1 763 764 ld_sstnight(nsstsets) = ln_sstnight 765 766 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 767 & sstfbfiles(jset:jset), & 768 & nsstvars, nsstextr, nitend-nit000+2, & 769 & dobsini, dobsend, ln_ignmis, .FALSE. ) 770 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 771 & ln_sst, ln_nea ) 772 773 END DO 774 775 ENDIF 776 777 ENDIF 778 779 ! - Sea surface salinity 780 IF ( ln_sss ) THEN 781 IF(lwp) WRITE(numout,*) ' SSS currently not available' 782 ENDIF 783 784 ! - Sea Ice Concentration 785 786 IF ( ln_seaice ) THEN 787 788 ! Set the number of variables for seaice to 1 789 nseaicevars = 1 790 791 ! Set the number of extra variables for seaice to 0 792 nseaiceextr = 0 793 794 ! Set the number of data sets to 1 795 nseaicesets = 1 796 797 ALLOCATE(seaicedata(nseaicesets)) 798 ALLOCATE(seaicedatqc(nseaicesets)) 799 seaicedata(:)%nsurf=0 800 seaicedatqc(:)%nsurf=0 801 802 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 803 & seaicefiles(1:jnumseaice), & 804 & nseaicevars, nseaiceextr, nitend-nit000+2, & 805 & dobsini, dobsend, ln_ignmis, .FALSE. ) 806 807 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 808 & ln_seaice, ln_nea ) 809 810 ENDIF 811 812 IF (ln_vel3d) THEN 813 814 ! Set the number of variables for profiles to 2 (U and V) 815 nvelovars = 2 816 817 ! Set the number of extra variables for profiles to 2 to store 818 ! rotation parameters 819 nveloextr = 2 820 821 jveloset = 0 822 823 IF ( ln_velavcur ) jveloset = jveloset + 1 824 IF ( ln_velhrcur ) jveloset = jveloset + 1 825 IF ( ln_velavadcp ) jveloset = jveloset + 1 826 IF ( ln_velhradcp ) jveloset = jveloset + 1 827 IF (ln_velfb) jveloset = jveloset + jnumvelfb 828 829 nvelosets = jveloset 830 IF ( nvelosets > 0 ) THEN 831 ALLOCATE( velodata(nvelosets) ) 832 ALLOCATE( veldatqc(nvelosets) ) 833 ALLOCATE( ld_velav(nvelosets) ) 834 ENDIF 835 836 jveloset = 0 837 838 ! Daily averaged data 839 840 IF ( ln_velavcur ) THEN 841 842 jveloset = jveloset + 1 843 844 ld_velav(jveloset) = .TRUE. 845 846 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 847 & velavcurfiles(1:jnumvelavcur), & 848 & nvelovars, nveloextr, & 849 & nitend-nit000+2, & 850 & dobsini, dobsend, ln_ignmis, & 851 & ld_velav(jveloset), & 852 & .FALSE. ) 853 854 DO jvar = 1, 2 855 CALL obs_prof_staend( velodata(jveloset), jvar ) 856 END DO 857 858 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 859 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 860 861 ENDIF 862 863 ! High frequency data 864 865 IF ( ln_velhrcur ) THEN 866 867 jveloset = jveloset + 1 868 869 ld_velav(jveloset) = .FALSE. 870 871 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 872 & velhrcurfiles(1:jnumvelhrcur), & 873 & nvelovars, nveloextr, & 874 & nitend-nit000+2, & 875 & dobsini, dobsend, ln_ignmis, & 876 & ld_velav(jveloset), & 877 & .FALSE. ) 878 879 DO jvar = 1, 2 880 CALL obs_prof_staend( velodata(jveloset), jvar ) 881 END DO 882 883 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 884 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 885 886 ENDIF 887 888 ! Daily averaged data 889 890 IF ( ln_velavadcp ) THEN 891 892 jveloset = jveloset + 1 893 894 ld_velav(jveloset) = .TRUE. 895 896 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 897 & velavadcpfiles(1:jnumvelavadcp), & 898 & nvelovars, nveloextr, & 899 & nitend-nit000+2, & 900 & dobsini, dobsend, ln_ignmis, & 901 & ld_velav(jveloset), & 902 & .FALSE. ) 903 904 DO jvar = 1, 2 905 CALL obs_prof_staend( velodata(jveloset), jvar ) 906 END DO 907 908 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 909 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 910 911 ENDIF 912 913 ! High frequency data 914 915 IF ( ln_velhradcp ) THEN 916 917 jveloset = jveloset + 1 918 919 ld_velav(jveloset) = .FALSE. 920 921 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 922 & velhradcpfiles(1:jnumvelhradcp), & 923 & nvelovars, nveloextr, & 924 & nitend-nit000+2, & 925 & dobsini, dobsend, ln_ignmis, & 926 & ld_velav(jveloset), & 927 & .FALSE. ) 928 929 DO jvar = 1, 2 930 CALL obs_prof_staend( velodata(jveloset), jvar ) 931 END DO 932 933 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 934 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 935 936 ENDIF 937 938 IF ( ln_velfb ) THEN 939 940 DO jset = 1, jnumvelfb 941 942 jveloset = jveloset + 1 943 944 ld_velav(jveloset) = ln_velfb_av(jset) 945 946 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 947 & velfbfiles(jset:jset), & 948 & nvelovars, nveloextr, & 949 & nitend-nit000+2, & 950 & dobsini, dobsend, ln_ignmis, & 951 & ld_velav(jveloset), & 952 & .FALSE. ) 953 954 DO jvar = 1, 2 955 CALL obs_prof_staend( velodata(jveloset), jvar ) 956 END DO 957 958 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 959 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 960 961 962 END DO 963 964 ENDIF 965 966 ENDIF 967 579 580 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 581 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 582 583 ENDIF 584 585 END DO 586 587 DEALLOCATE( ifilessurf, clsurffiles ) 588 589 ENDIF 590 591 CALL wrk_dealloc( jpi, jpj, zglam1 ) 592 CALL wrk_dealloc( jpi, jpj, zglam2 ) 593 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 594 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 595 CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 596 CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 597 968 598 END SUBROUTINE dia_obs_init 969 599 … … 975 605 !! 976 606 !! ** Method : Call the observation operators on each time step to 977 !! compute the model equivalent of the following date: 978 !! - T profiles 979 !! - S profiles 980 !! - Sea surface height (referenced to a mean) 981 !! - Sea surface temperature 982 !! - Sea surface salinity 983 !! - Velocity component (U,V) profiles 984 !! 985 !! ** Action : 607 !! compute the model equivalent of the following data: 608 !! - Profile data, currently T/S or U/V 609 !! - Surface data, currently SST, SLA or sea-ice concentration. 610 !! 611 !! ** Action : 986 612 !! 987 613 !! History : … … 992 618 !! ! 07-04 (G. Smith) Generalized surface operators 993 619 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 620 !! ! 15-08 (M. Martin) Combined surface/profile routines. 994 621 !!---------------------------------------------------------------------- 995 622 !! * Modules used 996 USE dom_oce, ONLY : & ! Ocean space and time domain variables 997 & rdt, & 998 & gdept_1d, & 999 & tmask, umask, vmask 1000 USE phycst, ONLY : & ! Physical constants 1001 & rday 1002 USE oce, ONLY : & ! Ocean dynamics and tracers variables 1003 & tsn, & 1004 & un, vn, & 623 USE phycst, ONLY : & ! Physical constants 624 & rday 625 USE oce, ONLY : & ! Ocean dynamics and tracers variables 626 & tsn, & 627 & un, & 628 & vn, & 1005 629 & sshn 1006 630 #if defined key_lim3 1007 USE ice, ONLY : & ! LIMIce model variables631 USE ice, ONLY : & ! LIM3 Ice model variables 1008 632 & frld 1009 633 #endif 1010 634 #if defined key_lim2 1011 USE ice_2, ONLY : & ! LIMIce model variables635 USE ice_2, ONLY : & ! LIM2 Ice model variables 1012 636 & frld 1013 637 #endif 638 #if defined key_cice 639 USE sbc_oce, ONLY : fr_i ! ice fraction 640 #endif 641 #if defined key_hadocc 642 USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2 643 & HADOCC_CHL, & 644 & HADOCC_FCO2, & 645 & HADOCC_PCO2, & 646 & HADOCC_FILL_FLT 647 #elif defined key_medusa && defined key_foam_medusa 648 USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2 649 & trn 650 USE par_medusa, ONLY: & 651 & jpchn, & 652 & jpchd 653 #if defined key_roam 654 USE sms_medusa, ONLY: & 655 & f2_pco2w, & 656 & f2_fco2w 657 #endif 658 #elif defined key_fabm 659 USE fabm 660 USE par_fabm 661 #endif 662 #if defined key_spm 663 USE par_spm, ONLY: & ! ERSEM/SPM sediments 664 & jp_spm 665 USE trc, ONLY : & 666 & trn 667 #endif 668 1014 669 IMPLICIT NONE 1015 670 1016 671 !! * Arguments 1017 INTEGER, INTENT(IN) :: kstp 672 INTEGER, INTENT(IN) :: kstp ! Current timestep 1018 673 !! * Local declarations 1019 INTEGER :: idaystp ! Number of timesteps per day 1020 INTEGER :: jprofset ! Profile data set loop variable 1021 INTEGER :: jslaset ! SLA data set loop variable 1022 INTEGER :: jsstset ! SST data set loop variable 1023 INTEGER :: jseaiceset ! sea ice data set loop variable 1024 INTEGER :: jveloset ! velocity profile data loop variable 1025 INTEGER :: jvar ! Variable number 1026 #if ! defined key_lim2 && ! defined key_lim3 1027 REAL(wp), POINTER, DIMENSION(:,:) :: frld 1028 #endif 1029 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1030 1031 #if ! defined key_lim2 && ! defined key_lim3 1032 CALL wrk_alloc(jpi,jpj,frld) 1033 #endif 674 INTEGER :: idaystp ! Number of timesteps per day 675 INTEGER :: jtype ! Data loop variable 676 INTEGER :: jvar ! Variable number 677 INTEGER :: ji, jj ! Loop counters 678 REAL(wp) :: tiny ! small number 679 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 680 & zprofvar1, & ! Model values for 1st variable in a prof ob 681 & zprofvar2 ! Model values for 2nd variable in a prof ob 682 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 683 & zprofmask1, & ! Mask associated with zprofvar1 684 & zprofmask2 ! Mask associated with zprofvar2 685 REAL(wp), POINTER, DIMENSION(:,:) :: & 686 & zsurfvar, & ! Model values equivalent to surface ob. 687 & zsurfmask ! Mask associated with surface variable 688 REAL(wp), POINTER, DIMENSION(:,:) :: & 689 & zglam1, & ! Model longitudes for prof variable 1 690 & zglam2, & ! Model longitudes for prof variable 2 691 & zgphi1, & ! Model latitudes for prof variable 1 692 & zgphi2 ! Model latitudes for prof variable 2 693 694 695 !Allocate local work arrays 696 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 697 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 698 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 699 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 700 CALL wrk_alloc( jpi, jpj, zsurfvar ) 701 CALL wrk_alloc( jpi, jpj, zsurfmask ) 702 CALL wrk_alloc( jpi, jpj, zglam1 ) 703 CALL wrk_alloc( jpi, jpj, zglam2 ) 704 CALL wrk_alloc( jpi, jpj, zgphi1 ) 705 CALL wrk_alloc( jpi, jpj, zgphi2 ) 1034 706 1035 707 IF(lwp) THEN … … 1037 709 WRITE(numout,*) 'dia_obs : Call the observation operators', kstp 1038 710 WRITE(numout,*) '~~~~~~~' 711 CALL FLUSH(numout) 1039 712 ENDIF 1040 713 … … 1042 715 1043 716 !----------------------------------------------------------------------- 1044 ! No LIM => frld == 0.0_wp717 ! Call the profile and surface observation operators 1045 718 !----------------------------------------------------------------------- 1046 #if ! defined key_lim2 && ! defined key_lim3 1047 frld(:,:) = 0.0_wp 719 720 IF ( nproftypes > 0 ) THEN 721 722 DO jtype = 1, nproftypes 723 724 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 725 CASE('prof') 726 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 727 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 728 zprofmask1(:,:,:) = tmask(:,:,:) 729 zprofmask2(:,:,:) = tmask(:,:,:) 730 zglam1(:,:) = glamt(:,:) 731 zglam2(:,:) = glamt(:,:) 732 zgphi1(:,:) = gphit(:,:) 733 zgphi2(:,:) = gphit(:,:) 734 CASE('vel') 735 zprofvar1(:,:,:) = un(:,:,:) 736 zprofvar2(:,:,:) = vn(:,:,:) 737 zprofmask1(:,:,:) = umask(:,:,:) 738 zprofmask2(:,:,:) = vmask(:,:,:) 739 zglam1(:,:) = glamu(:,:) 740 zglam2(:,:) = glamv(:,:) 741 zgphi1(:,:) = gphiu(:,:) 742 zgphi2(:,:) = gphiv(:,:) 743 CASE DEFAULT 744 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 745 END SELECT 746 747 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 748 & nit000, idaystp, & 749 & zprofvar1, zprofvar2, & 750 & fsdept(:,:,:), fsdepw(:,:,:), & 751 & zprofmask1, zprofmask2, & 752 & zglam1, zglam2, zgphi1, zgphi2, & 753 & nn_1dint, nn_2dint, & 754 & kdailyavtypes = nn_profdavtypes ) 755 756 END DO 757 758 ENDIF 759 760 IF ( nsurftypes > 0 ) THEN 761 762 DO jtype = 1, nsurftypes 763 764 !Defaults which might be changed 765 zsurfmask(:,:) = tmask(:,:,1) 766 767 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 768 CASE('sst') 769 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 770 CASE('sla') 771 zsurfvar(:,:) = sshn(:,:) 772 CASE('sss') 773 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 774 CASE('sic') 775 IF ( kstp == 0 ) THEN 776 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 777 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 778 & 'time-step but some obs are valid then.' ) 779 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 780 & ' sea-ice obs will be missed' 781 ENDIF 782 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 783 & surfdataqc(jtype)%nsstp(1) 784 CYCLE 785 ELSE 786 #if defined key_cice 787 zsurfvar(:,:) = fr_i(:,:) 788 #elif defined key_lim2 || defined key_lim3 789 zsurfvar(:,:) = 1._wp - frld(:,:) 790 #else 791 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 792 & ' but no sea-ice model appears to have been defined' ) 1048 793 #endif 1049 !----------------------------------------------------------------------- 1050 ! Depending on switches call various observation operators 1051 !----------------------------------------------------------------------- 1052 1053 ! - Temperature/salinity profiles 1054 IF ( ln_t3d .OR. ln_s3d ) THEN 1055 DO jprofset = 1, nprofsets 1056 IF ( ld_enact(jprofset) ) THEN 1057 CALL obs_pro_opt( prodatqc(jprofset), & 1058 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1059 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1060 & gdept_1d, tmask, n1dint, n2dint, & 1061 & kdailyavtypes = endailyavtypes ) 1062 ELSE 1063 CALL obs_pro_opt( prodatqc(jprofset), & 1064 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1065 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1066 & gdept_1d, tmask, n1dint, n2dint ) 1067 ENDIF 794 ENDIF 795 796 CASE('logchl') 797 #if defined key_hadocc 798 zsurfvar(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC 799 #elif defined key_medusa && defined key_foam_medusa 800 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 801 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 802 #elif defined key_fabm 803 chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot) 804 zsurfvar(:,:) = chl_3d(:,:,1) 805 #else 806 CALL ctl_stop( ' Trying to run logchl observation operator', & 807 & ' but no biogeochemical model appears to have been defined' ) 808 #endif 809 zsurfmask(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things 810 ! Take the log10 where we can, otherwise exclude 811 tiny = 1.0e-20 812 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 813 zsurfvar(:,:) = LOG10(zsurfvar(:,:)) 814 ELSEWHERE 815 zsurfvar(:,:) = obfillflt 816 zsurfmask(:,:) = 0 817 END WHERE 818 CASE('spm') 819 #if defined key_spm 820 zsurfvar(:,:) = 0.0 821 DO jn = 1, jp_spm 822 zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes 823 END DO 824 #else 825 CALL ctl_stop( ' Trying to run spm observation operator', & 826 & ' but no spm model appears to have been defined' ) 827 #endif 828 CASE('fco2') 829 #if defined key_hadocc 830 zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 831 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 832 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 833 zsurfvar(:,:) = obfillflt 834 zsurfmask(:,:) = 0 835 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 836 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 837 ENDIF 838 #elif defined key_medusa && defined key_foam_medusa && defined key_roam 839 zsurfvar(:,:) = f2_fco2w(:,:) 840 #elif defined key_fabm 841 ! First, get pCO2 from FABM 842 pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 843 zsurfvar(:,:) = pco2_3d(:,:,1) 844 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 845 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 846 ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 847 ! and 848 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 849 ! Marine Chemistry, 2: 203-215. 850 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 851 ! not explicitly included - atmospheric pressure is not necessarily available so this is 852 ! the best assumption. 853 ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 854 ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 855 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 856 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 857 zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & 858 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & 859 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 860 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 861 & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 862 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 863 #else 864 CALL ctl_stop( ' Trying to run fco2 observation operator', & 865 & ' but no biogeochemical model appears to have been defined' ) 866 #endif 867 CASE('pco2') 868 #if defined key_hadocc 869 zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 870 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 871 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 872 zsurfvar(:,:) = obfillflt 873 zsurfmask(:,:) = 0 874 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 875 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 876 ENDIF 877 #elif defined key_medusa && defined key_foam_medusa && defined key_roam 878 zsurfvar(:,:) = f2_pco2w(:,:) 879 #elif defined key_fabm 880 pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 881 zsurfvar(:,:) = pco2_3d(:,:,1) 882 #else 883 CALL ctl_stop( ' Trying to run pCO2 observation operator', & 884 & ' but no biogeochemical model appears to have been defined' ) 885 #endif 886 887 CASE DEFAULT 888 889 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 890 891 END SELECT 892 893 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 894 & nit000, idaystp, zsurfvar, zsurfmask, & 895 & n2dintsurf(jtype), llnightav(jtype), & 896 & ravglamscl(jtype), ravgphiscl(jtype), & 897 & lfpindegs(jtype) ) 898 1068 899 END DO 1069 ENDIF 1070 1071 ! - Sea surface anomaly 1072 IF ( ln_sla ) THEN 1073 DO jslaset = 1, nslasets 1074 CALL obs_sla_opt( sladatqc(jslaset), & 1075 & kstp, jpi, jpj, nit000, sshn, & 1076 & tmask(:,:,1), n2dint ) 1077 END DO 1078 ENDIF 1079 1080 ! - Sea surface temperature 1081 IF ( ln_sst ) THEN 1082 DO jsstset = 1, nsstsets 1083 CALL obs_sst_opt( sstdatqc(jsstset), & 1084 & kstp, jpi, jpj, nit000, idaystp, & 1085 & tsn(:,:,1,jp_tem), tmask(:,:,1), & 1086 & n2dint, ld_sstnight(jsstset) ) 1087 END DO 1088 ENDIF 1089 1090 ! - Sea surface salinity 1091 IF ( ln_sss ) THEN 1092 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1093 ENDIF 1094 1095 #if defined key_lim2 || defined key_lim3 1096 IF ( ln_seaice ) THEN 1097 DO jseaiceset = 1, nseaicesets 1098 CALL obs_seaice_opt( seaicedatqc(jseaiceset), & 1099 & kstp, jpi, jpj, nit000, 1.-frld, & 1100 & tmask(:,:,1), n2dint ) 1101 END DO 1102 ENDIF 1103 #endif 1104 1105 ! - Velocity profiles 1106 IF ( ln_vel3d ) THEN 1107 DO jveloset = 1, nvelosets 1108 ! zonal component of velocity 1109 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, & 1110 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, & 1111 n1dint, n2dint, ld_velav(jveloset) ) 1112 END DO 1113 ENDIF 1114 1115 #if ! defined key_lim2 && ! defined key_lim3 1116 CALL wrk_dealloc(jpi,jpj,frld) 1117 #endif 900 901 ENDIF 902 903 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 904 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 905 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 906 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 907 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 908 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 909 CALL wrk_dealloc( jpi, jpj, zglam1 ) 910 CALL wrk_dealloc( jpi, jpj, zglam2 ) 911 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 912 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 1118 913 1119 914 END SUBROUTINE dia_obs 1120 1121 SUBROUTINE dia_obs_wri 915 916 SUBROUTINE dia_obs_wri 1122 917 !!---------------------------------------------------------------------- 1123 918 !! *** ROUTINE dia_obs_wri *** … … 1127 922 !! ** Method : Call observation diagnostic output routines 1128 923 !! 1129 !! ** Action : 924 !! ** Action : 1130 925 !! 1131 926 !! History : … … 1135 930 !! ! 07-03 (K. Mogensen) General handling of profiles 1136 931 !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles 1137 !!---------------------------------------------------------------------- 932 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 933 !!---------------------------------------------------------------------- 934 !! * Modules used 935 USE obs_rot_vel ! Rotation of velocities 936 1138 937 IMPLICIT NONE 1139 938 1140 939 !! * Local declarations 1141 1142 INTEGER :: jprofset ! Profile data set loop variable 1143 INTEGER :: jveloset ! Velocity data set loop variable 1144 INTEGER :: jslaset ! SLA data set loop variable 1145 INTEGER :: jsstset ! SST data set loop variable 1146 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1147 INTEGER :: jset 1148 INTEGER :: jfbini 1149 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1150 CHARACTER(LEN=10) :: cdtmp 940 INTEGER :: jtype ! Data set loop variable 941 INTEGER :: jo, jvar, jk 942 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 943 & zu, & 944 & zv 945 1151 946 !----------------------------------------------------------------------- 1152 947 ! Depending on switches call various observation output routines 1153 948 !----------------------------------------------------------------------- 1154 949 1155 ! - Temperature/salinity profiles 1156 1157 IF( ln_t3d .OR. ln_s3d ) THEN 1158 1159 ! Copy data from prodatqc to profdata structures 1160 DO jprofset = 1, nprofsets 1161 1162 CALL obs_prof_decompress( prodatqc(jprofset), & 1163 & profdata(jprofset), .TRUE., numout ) 950 IF ( nproftypes > 0 ) THEN 951 952 DO jtype = 1, nproftypes 953 954 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 955 956 ! For velocity data, rotate the model velocities to N/S, E/W 957 ! using the compressed data structure. 958 ALLOCATE( & 959 & zu(profdataqc(jtype)%nvprot(1)), & 960 & zv(profdataqc(jtype)%nvprot(2)) & 961 & ) 962 963 CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 964 965 DO jo = 1, profdataqc(jtype)%nprof 966 DO jvar = 1, 2 967 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 968 969 IF ( jvar == 1 ) THEN 970 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 971 ELSE 972 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 973 ENDIF 974 975 END DO 976 END DO 977 END DO 978 979 DEALLOCATE( zu ) 980 DEALLOCATE( zv ) 981 982 END IF 983 984 CALL obs_prof_decompress( profdataqc(jtype), & 985 & profdata(jtype), .TRUE., numout ) 986 987 CALL obs_wri_prof( profdata(jtype) ) 1164 988 1165 989 END DO 1166 990 1167 ! Write the profiles. 1168 1169 jprofset = 0 1170 1171 ! ENACT insitu data 1172 1173 IF ( ln_ena ) THEN 1174 1175 jprofset = jprofset + 1 1176 1177 CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 1178 1179 ENDIF 1180 1181 ! Coriolis insitu data 1182 1183 IF ( ln_cor ) THEN 1184 1185 jprofset = jprofset + 1 1186 1187 CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 1188 1189 ENDIF 1190 1191 ! Feedback insitu data 1192 1193 IF ( ln_profb ) THEN 1194 1195 jfbini = jprofset + 1 1196 1197 DO jprofset = jfbini, nprofsets 1198 1199 jset = jprofset - jfbini + 1 1200 WRITE(cdtmp,'(A,I2.2)')'profb_',jset 1201 CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 1202 1203 END DO 1204 1205 ENDIF 1206 1207 ENDIF 1208 1209 ! - Sea surface anomaly 1210 IF ( ln_sla ) THEN 1211 1212 ! Copy data from sladatqc to sladata structures 1213 DO jslaset = 1, nslasets 1214 1215 CALL obs_surf_decompress( sladatqc(jslaset), & 1216 & sladata(jslaset), .TRUE., numout ) 991 ENDIF 992 993 IF ( nsurftypes > 0 ) THEN 994 995 DO jtype = 1, nsurftypes 996 997 CALL obs_surf_decompress( surfdataqc(jtype), & 998 & surfdata(jtype), .TRUE., numout ) 999 1000 CALL obs_wri_surf( surfdata(jtype) ) 1217 1001 1218 1002 END DO 1219 1003 1220 jslaset = 01221 1222 ! Write the AVISO SLA data1223 1224 IF ( ln_sladt ) THEN1225 1226 jslaset = 11227 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )1228 jslaset = 21229 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )1230 1231 ENDIF1232 1233 IF ( ln_slafb ) THEN1234 1235 jfbini = jslaset + 11236 1237 DO jslaset = jfbini, nslasets1238 1239 jset = jslaset - jfbini + 11240 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset1241 CALL obs_wri_sla( cdtmp, sladata(jslaset) )1242 1243 END DO1244 1245 ENDIF1246 1247 ENDIF1248 1249 ! - Sea surface temperature1250 IF ( ln_sst ) THEN1251 1252 ! Copy data from sstdatqc to sstdata structures1253 DO jsstset = 1, nsstsets1254 1255 CALL obs_surf_decompress( sstdatqc(jsstset), &1256 & sstdata(jsstset), .TRUE., numout )1257 1258 END DO1259 1260 jsstset = 01261 1262 ! Write the AVISO SST data1263 1264 IF ( ln_reysst ) THEN1265 1266 jsstset = jsstset + 11267 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )1268 1269 ENDIF1270 1271 IF ( ln_ghrsst ) THEN1272 1273 jsstset = jsstset + 11274 CALL obs_wri_sst( 'ghr', sstdata(jsstset) )1275 1276 ENDIF1277 1278 IF ( ln_sstfb ) THEN1279 1280 jfbini = jsstset + 11281 1282 DO jsstset = jfbini, nsstsets1283 1284 jset = jsstset - jfbini + 11285 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset1286 CALL obs_wri_sst( cdtmp, sstdata(jsstset) )1287 1288 END DO1289 1290 ENDIF1291 1292 ENDIF1293 1294 ! - Sea surface salinity1295 IF ( ln_sss ) THEN1296 IF(lwp) WRITE(numout,*) ' SSS currently not available'1297 ENDIF1298 1299 ! - Sea Ice Concentration1300 IF ( ln_seaice ) THEN1301 1302 ! Copy data from seaicedatqc to seaicedata structures1303 DO jseaiceset = 1, nseaicesets1304 1305 CALL obs_surf_decompress( seaicedatqc(jseaiceset), &1306 & seaicedata(jseaiceset), .TRUE., numout )1307 1308 END DO1309 1310 ! Write the Sea Ice data1311 DO jseaiceset = 1, nseaicesets1312 1313 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset1314 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )1315 1316 END DO1317 1318 ENDIF1319 1320 ! Velocity data1321 IF( ln_vel3d ) THEN1322 1323 ! Copy data from veldatqc to velodata structures1324 DO jveloset = 1, nvelosets1325 1326 CALL obs_prof_decompress( veldatqc(jveloset), &1327 & velodata(jveloset), .TRUE., numout )1328 1329 END DO1330 1331 ! Write the profiles.1332 1333 jveloset = 01334 1335 ! Daily averaged data1336 1337 IF ( ln_velavcur ) THEN1338 1339 jveloset = jveloset + 11340 1341 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )1342 1343 ENDIF1344 1345 ! High frequency data1346 1347 IF ( ln_velhrcur ) THEN1348 1349 jveloset = jveloset + 11350 1351 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )1352 1353 ENDIF1354 1355 ! Daily averaged data1356 1357 IF ( ln_velavadcp ) THEN1358 1359 jveloset = jveloset + 11360 1361 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )1362 1363 ENDIF1364 1365 ! High frequency data1366 1367 IF ( ln_velhradcp ) THEN1368 1369 jveloset = jveloset + 11370 1371 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )1372 1373 ENDIF1374 1375 ! Feedback velocity data1376 1377 IF ( ln_velfb ) THEN1378 1379 jfbini = jveloset + 11380 1381 DO jveloset = jfbini, nvelosets1382 1383 jset = jveloset - jfbini + 11384 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset1385 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )1386 1387 END DO1388 1389 ENDIF1390 1391 1004 ENDIF 1392 1005 … … 1406 1019 !! 1407 1020 !!---------------------------------------------------------------------- 1408 ! !obs_grid deallocation1021 ! obs_grid deallocation 1409 1022 CALL obs_grid_deallocate 1410 1023 1411 !! diaobs deallocation 1412 IF ( nprofsets > 0 ) THEN 1413 DEALLOCATE(ld_enact, & 1414 & profdata, & 1415 & prodatqc) 1416 END IF 1417 IF ( ln_sla ) THEN 1418 DEALLOCATE(sladata, & 1419 & sladatqc) 1420 END IF 1421 IF ( ln_seaice ) THEN 1422 DEALLOCATE(sladata, & 1423 & sladatqc) 1424 END IF 1425 IF ( ln_sst ) THEN 1426 DEALLOCATE(sstdata, & 1427 & sstdatqc) 1428 END IF 1429 IF ( ln_vel3d ) THEN 1430 DEALLOCATE(ld_velav, & 1431 & velodata, & 1432 & veldatqc) 1433 END IF 1024 ! diaobs deallocation 1025 IF ( nproftypes > 0 ) & 1026 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 1027 1028 IF ( nsurftypes > 0 ) & 1029 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 1030 & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 1031 1434 1032 END SUBROUTINE dia_obs_dealloc 1435 1033 … … 1454 1052 USE phycst, ONLY : & ! Physical constants 1455 1053 & rday 1456 ! USE daymod, ONLY : & ! Time variables1457 ! & nmonth_len1458 1054 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1459 1055 & rdt … … 1472 1068 INTEGER :: imin 1473 1069 INTEGER :: imday ! Number of days in month. 1474 REAL( KIND=wp) :: zdayfrc ! Fraction of day1070 REAL(wp) :: zdayfrc ! Fraction of day 1475 1071 1476 1072 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year 1477 1073 1478 ! !----------------------------------------------------------------------1479 ! !Initial date initialization (year, month, day, hour, minute)1480 ! !----------------------------------------------------------------------1074 !---------------------------------------------------------------------- 1075 ! Initial date initialization (year, month, day, hour, minute) 1076 !---------------------------------------------------------------------- 1481 1077 iyea = ndate0 / 10000 1482 1078 imon = ( ndate0 - iyea * 10000 ) / 100 … … 1485 1081 imin = ( nn_time0 - ihou * 100 ) 1486 1082 1487 ! !----------------------------------------------------------------------1488 ! !Compute number of days + number of hours + min since initial time1489 ! !----------------------------------------------------------------------1083 !---------------------------------------------------------------------- 1084 ! Compute number of days + number of hours + min since initial time 1085 !---------------------------------------------------------------------- 1490 1086 zdayfrc = kstp * rdt / rday 1491 1087 zdayfrc = zdayfrc - aint(zdayfrc) … … 1501 1097 iday = iday + kstp * rdt / rday 1502 1098 1503 ! !-----------------------------------------------------------------------1504 ! !Convert number of days (iday) into a real date1505 ! !----------------------------------------------------------------------1099 !----------------------------------------------------------------------- 1100 ! Convert number of days (iday) into a real date 1101 !---------------------------------------------------------------------- 1506 1102 1507 1103 CALL calc_month_len( iyea, imonth_len ) … … 1517 1113 END DO 1518 1114 1519 ! !----------------------------------------------------------------------1520 ! !Convert it into YYYYMMDD.HHMMSS format.1521 ! !----------------------------------------------------------------------1115 !---------------------------------------------------------------------- 1116 ! Convert it into YYYYMMDD.HHMMSS format. 1117 !---------------------------------------------------------------------- 1522 1118 ddobs = iyea * 10000_dp + imon * 100_dp + & 1523 1119 & iday + ihou * 0.01_dp + imin * 0.0001_dp … … 1528 1124 !!---------------------------------------------------------------------- 1529 1125 !! *** ROUTINE ini_date *** 1530 !! 1531 !! ** Purpose : Get initial dat ain double precision YYYYMMDD.HHMMSS format1532 !! 1533 !! ** Method : Get initial dat ain double precision YYYYMMDD.HHMMSS format1534 !! 1535 !! ** Action : Get initial dat ain double precision YYYYMMDD.HHMMSS format1126 !! 1127 !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 1128 !! 1129 !! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format 1130 !! 1131 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1536 1132 !! 1537 1133 !! History : … … 1546 1142 1547 1143 !! * Arguments 1548 REAL( KIND=dp), INTENT(OUT) :: ddobsini! Initial date in YYYYMMDD.HHMMSS1144 REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 1549 1145 1550 1146 CALL calc_date( nit000 - 1, ddobsini ) … … 1555 1151 !!---------------------------------------------------------------------- 1556 1152 !! *** ROUTINE fin_date *** 1557 !! 1558 !! ** Purpose : Get final dat ain double precision YYYYMMDD.HHMMSS format1559 !! 1560 !! ** Method : Get final dat ain double precision YYYYMMDD.HHMMSS format1561 !! 1562 !! ** Action : Get final dat ain double precision YYYYMMDD.HHMMSS format1153 !! 1154 !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 1155 !! 1156 !! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format 1157 !! 1158 !! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format 1563 1159 !! 1564 1160 !! History : … … 1572 1168 1573 1169 !! * Arguments 1574 REAL( KIND=dp), INTENT(OUT) :: ddobsfin! Final date in YYYYMMDD.HHMMSS1170 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1575 1171 1576 1172 CALL calc_date( nitend, ddobsfin ) 1577 1173 1578 END SUBROUTINE fin_date 1579 1174 END SUBROUTINE fin_date 1175 1176 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 1177 & cfilestype, ifiles, cobstypes, cfiles ) 1178 1179 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1180 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 1181 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1182 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1183 & ifiles ! Out appended number of files for this type 1184 1185 CHARACTER(len=6), INTENT(IN) :: ctypein 1186 CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 1187 & cfilestype ! In list of files for this obs type 1188 CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 1189 & cobstypes ! Out appended list of obs types 1190 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 1191 & cfiles ! Out appended list of files for all types 1192 1193 !Local variables 1194 INTEGER :: jfile 1195 1196 cfiles(jtype,:) = cfilestype(:) 1197 cobstypes(jtype) = ctypein 1198 ifiles(jtype) = 0 1199 DO jfile = 1, jpmaxnfiles 1200 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 1201 ifiles(jtype) = ifiles(jtype) + 1 1202 END DO 1203 1204 IF ( ifiles(jtype) == 0 ) THEN 1205 CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// & 1206 & ' set to true but no files available to read' ) 1207 ENDIF 1208 1209 IF(lwp) THEN 1210 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1211 DO jfile = 1, ifiles(jtype) 1212 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1213 END DO 1214 ENDIF 1215 1216 END SUBROUTINE obs_settypefiles 1217 1218 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1219 & n2dint_default, n2dint_type, & 1220 & ravglamscl_type, ravgphiscl_type, & 1221 & lfp_indegs_type, lavnight_type, & 1222 & n2dint, ravglamscl, ravgphiscl, & 1223 & lfpindegs, lavnight ) 1224 1225 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1226 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1227 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1228 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1229 REAL(wp), INTENT(IN) :: & 1230 & ravglamscl_type, & !E/W diameter of obs footprint for this type 1231 & ravgphiscl_type !N/S diameter of obs footprint for this type 1232 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1233 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1234 CHARACTER(len=6), INTENT(IN) :: ctypein 1235 1236 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1237 & n2dint 1238 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 1239 & ravglamscl, ravgphiscl 1240 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 1241 & lfpindegs, lavnight 1242 1243 lavnight(jtype) = lavnight_type 1244 1245 IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 1246 n2dint(jtype) = n2dint_type 1247 ELSE 1248 n2dint(jtype) = n2dint_default 1249 ENDIF 1250 1251 ! For averaging observation footprints set options for size of footprint 1252 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 1253 IF ( ravglamscl_type > 0._wp ) THEN 1254 ravglamscl(jtype) = ravglamscl_type 1255 ELSE 1256 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1257 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) 1258 ENDIF 1259 1260 IF ( ravgphiscl_type > 0._wp ) THEN 1261 ravgphiscl(jtype) = ravgphiscl_type 1262 ELSE 1263 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1264 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) 1265 ENDIF 1266 1267 lfpindegs(jtype) = lfp_indegs_type 1268 1269 ENDIF 1270 1271 ! Write out info 1272 IF(lwp) THEN 1273 IF ( n2dint(jtype) <= 4 ) THEN 1274 WRITE(numout,*) ' '//TRIM(ctypein)// & 1275 & ' model counterparts will be interpolated horizontally' 1276 ELSE IF ( n2dint(jtype) <= 6 ) THEN 1277 WRITE(numout,*) ' '//TRIM(ctypein)// & 1278 & ' model counterparts will be averaged horizontally' 1279 WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) 1280 WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) 1281 IF ( lfpindegs(jtype) ) THEN 1282 WRITE(numout,*) ' '//' (in degrees)' 1283 ELSE 1284 WRITE(numout,*) ' '//' (in metres)' 1285 ENDIF 1286 ENDIF 1287 ENDIF 1288 1289 END SUBROUTINE obs_setinterpopts 1290 1580 1291 END MODULE diaobs -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90
r2358 r10120 325 325 CALL obs_mpp_max_integer( kobsj, kobs ) 326 326 ELSE 327 CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj,kobs )327 CALL obs_mpp_find_obs_proc( kproc,kobs ) 328 328 ENDIF 329 329 -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90
r6486 r10120 52 52 53 53 !! Default values 54 REAL, PUBLIC :: grid_search_res = 0.5! Resolution of grid54 REAL, PUBLIC :: rn_gridsearchres = 0.5 ! Resolution of grid 55 55 INTEGER, PRIVATE :: gsearch_nlons_def ! Num of longitudes 56 56 INTEGER, PRIVATE :: gsearch_nlats_def ! Num of latitudes … … 83 83 LOGICAL, PUBLIC :: ln_grid_global ! Use global distribution of observations 84 84 CHARACTER(LEN=44), PUBLIC :: & 85 & grid_search_file ! file name head for grid search lookup85 & cn_gridsearchfile ! file name head for grid search lookup 86 86 87 87 !!---------------------------------------------------------------------- … … 613 613 CALL obs_mpp_max_integer( kobsj, kobs ) 614 614 ELSE 615 CALL obs_mpp_find_obs_proc( kproc, kobs i, kobsj, kobs)615 CALL obs_mpp_find_obs_proc( kproc, kobs ) 616 616 ENDIF 617 617 … … 690 690 691 691 IF(lwp) WRITE(numout,*) 692 IF(lwp) WRITE(numout,*)'Grid search resolution : ', grid_search_res693 694 gsearch_nlons_def = NINT( 360.0_wp / grid_search_res )695 gsearch_nlats_def = NINT( 180.0_wp / grid_search_res )696 gsearch_lonmin_def = -180.0_wp + 0.5_wp * grid_search_res697 gsearch_latmin_def = -90.0_wp + 0.5_wp * grid_search_res698 gsearch_dlon_def = grid_search_res699 gsearch_dlat_def = grid_search_res692 IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres 693 694 gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres ) 695 gsearch_nlats_def = NINT( 180.0_wp / rn_gridsearchres ) 696 gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres 697 gsearch_latmin_def = -90.0_wp + 0.5_wp * rn_gridsearchres 698 gsearch_dlon_def = rn_gridsearchres 699 gsearch_dlat_def = rn_gridsearchres 700 700 701 701 IF (lwp) THEN … … 710 710 IF ( ln_grid_global ) THEN 711 711 WRITE(cfname, FMT="(A,'_',A)") & 712 & TRIM( grid_search_file), 'global.nc'712 & TRIM(cn_gridsearchfile), 'global.nc' 713 713 ELSE 714 714 WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") & 715 & TRIM( grid_search_file), nproc, jpni, jpnj715 & TRIM(cn_gridsearchfile), nproc, jpni, jpnj 716 716 ENDIF 717 717 -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90
r6486 r10120 35 35 CONTAINS 36 36 37 SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &37 SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 38 38 & pval, pgval, kproc ) 39 39 !!---------------------------------------------------------------------- … … 57 57 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 58 58 INTEGER, INTENT(IN) :: kobs ! Local number of observations 59 INTEGER, INTENT(IN) :: kpi ! Number of points in i direction 60 INTEGER, INTENT(IN) :: kpj ! Number of points in j direction 59 61 INTEGER, INTENT(IN) :: kpk ! Number of levels 60 62 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & … … 63 65 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 64 66 & kproc ! Precomputed processor for each i,j,iobs points 65 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&67 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 66 68 & pval ! Local 3D array to extract data from 67 69 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& … … 73 75 IF (PRESENT(kproc)) THEN 74 76 75 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, &77 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 76 78 & kgrdj, pval, pgval, kproc=kproc ) 77 79 78 80 ELSE 79 81 80 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, &82 CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 81 83 & kgrdj, pval, pgval ) 82 84 … … 85 87 ELSE 86 88 87 CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &89 CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 88 90 & pval, pgval ) 89 91 … … 92 94 END SUBROUTINE obs_int_comm_3d 93 95 94 SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, k grdi, kgrdj, pval, pgval, &96 SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kpi, kpj, kgrdi, kgrdj, pval, pgval, & 95 97 & kproc ) 96 98 !!---------------------------------------------------------------------- … … 111 113 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 112 114 INTEGER, INTENT(IN) :: kobs ! Local number of observations 115 INTEGER, INTENT(IN) :: kpi ! Number of model grid points in i direction 116 INTEGER, INTENT(IN) :: kpj ! Number of model grid points in j direction 113 117 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 114 118 & kgrdi, & ! i,j indicies for each stencil … … 116 120 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 117 121 & kproc ! Precomputed processor for each i,j,iobs points 118 REAL(KIND=wp), DIMENSION( jpi,jpj), INTENT(IN) ::&122 REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) ::& 119 123 & pval ! Local 3D array to extra data from 120 124 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::& … … 136 140 IF (PRESENT(kproc)) THEN 137 141 138 CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &142 CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 139 143 & zgval, kproc=kproc ) 140 144 ELSE 141 145 142 CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &146 CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 143 147 & zgval ) 144 148 … … 154 158 END SUBROUTINE obs_int_comm_2d 155 159 156 SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &160 SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 157 161 & pval, pgval, kproc ) 158 162 !!---------------------------------------------------------------------- … … 174 178 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 175 179 INTEGER, INTENT(IN) :: kobs ! Local number of observations 180 INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction 181 INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction 176 182 INTEGER, INTENT(IN) :: kpk ! Number of levels 177 183 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & … … 180 186 INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 181 187 & kproc ! Precomputed processor for each i,j,iobs points 182 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&188 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 183 189 & pval ! Local 3D array to extract data from 184 190 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& … … 207 213 208 214 ! Check valid points 209 215 210 216 IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. & 211 217 & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN 212 218 213 219 CALL ctl_stop( 'Error in obs_int_comm_3d_global', & 214 220 & 'Point outside global domain' ) 215 221 216 222 ENDIF 217 223 … … 323 329 END SUBROUTINE obs_int_comm_3d_global 324 330 325 SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kp k, kgrdi, kgrdj, &331 SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 326 332 & pval, pgval ) 327 333 !!---------------------------------------------------------------------- … … 343 349 INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil 344 350 INTEGER, INTENT(IN) :: kobs ! Local number of observations 351 INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction 352 INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction 345 353 INTEGER, INTENT(IN) :: kpk ! Number of levels 346 354 INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 347 355 & kgrdi, & ! i,j indicies for each stencil 348 356 & kgrdj 349 REAL(KIND=wp), DIMENSION( jpi,jpj,kpk), INTENT(IN) ::&357 REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 350 358 & pval ! Local 3D array to extract data from 351 359 REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90
r6486 r10120 7 7 !! - ! 2006-05 (K. Mogensen) Reformatted 8 8 !! - ! 2008-01 (K. Mogensen) add mpp_global_max 9 !! 3.6 ! 2015-01 (J. Waters) obs_mpp_find_obs_proc 10 !! rewritten to avoid global arrays 9 11 !!---------------------------------------------------------------------- 10 12 # define mpivar mpi_double_precision … … 12 14 !! obs_mpp_bcast_integer : Broadcast an integer array from a processor to all processors 13 15 !! obs_mpp_max_integer : Find maximum on all processors of each value in an integer on all processors 14 !! obs_mpp_find_obs_proc : Find processors which should hold the observations 16 !! obs_mpp_find_obs_proc : Find processors which should hold the observations, avoiding global arrays 15 17 !! obs_mpp_sum_integers : Sum an integer array from all processors 16 18 !! obs_mpp_sum_integer : Sum an integer from all processors … … 96 98 ! 97 99 INTEGER :: ierr 98 INTEGER, DIMENSION(kno) :: ivals 99 ! 100 INCLUDE 'mpif.h' 101 !!---------------------------------------------------------------------- 100 INTEGER, DIMENSION(:), ALLOCATABLE :: ivals 101 ! 102 INCLUDE 'mpif.h' 103 !!---------------------------------------------------------------------- 104 105 ALLOCATE( ivals(kno) ) 102 106 103 107 ! Call the MPI library to find the maximum across processors … … 105 109 & mpi_max, mpi_comm_opa, ierr ) 106 110 kvals(:) = ivals(:) 111 112 DEALLOCATE( ivals ) 107 113 #else 108 114 ! no MPI: empty routine … … 111 117 112 118 113 SUBROUTINE obs_mpp_find_obs_proc( kobsp, kobsi, kobsj,kno )114 !!---------------------------------------------------------------------- 115 !! *** ROUTINE obs_mpp_find_obs_proc ***116 !! 117 !! ** Purpose : From the array kobsp containing the results of the grid119 SUBROUTINE obs_mpp_find_obs_proc( kobsp,kno ) 120 !!---------------------------------------------------------------------- 121 !! *** ROUTINE obs_mpp_find_obs_proc *** 122 !! 123 !! ** Purpose : From the array kobsp containing the results of the 118 124 !! grid search on each processor the processor return a 119 125 !! decision of which processors should hold the observation. 120 126 !! 121 !! ** Method : A temporary 2D array holding all the decisions is122 !! constructed using mpi_allgather on each processor.123 !! If more than one processor has found the observation124 !! with the observation in the inner domain gets it125 !! 126 !! ** Action : This does only work for MPI. 127 !! ** Method : Synchronize the processor number for each obs using 128 !! obs_mpp_max_integer. If an observation exists on two 129 !! processors it will be allocated to the lower numbered 130 !! processor. 131 !! 132 !! ** Action : This does only work for MPI. 127 133 !! It does not work for SHMEM. 128 134 !! … … 130 136 !!---------------------------------------------------------------------- 131 137 INTEGER , INTENT(in ) :: kno 132 INTEGER, DIMENSION(kno), INTENT(in ) :: kobsi, kobsj133 138 INTEGER, DIMENSION(kno), INTENT(inout) :: kobsp 134 139 ! 135 140 #if defined key_mpp_mpi 136 141 ! 137 INTEGER :: ji 138 INTEGER :: jj 139 INTEGER :: size 140 INTEGER :: ierr 141 INTEGER :: iobsip 142 INTEGER :: iobsjp 143 INTEGER :: num_sus_obs 144 INTEGER, DIMENSION(kno) :: iobsig, iobsjg 145 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iobsp, iobsi, iobsj 146 !! 147 INCLUDE 'mpif.h' 148 !!---------------------------------------------------------------------- 149 150 !----------------------------------------------------------------------- 151 ! Call the MPI library to find the maximum accross processors 152 !----------------------------------------------------------------------- 153 CALL mpi_comm_size( mpi_comm_opa, size, ierr ) 154 !----------------------------------------------------------------------- 155 ! Convert local grids points to global grid points 156 !----------------------------------------------------------------------- 142 ! 143 INTEGER :: ji, isum 144 INTEGER, DIMENSION(:), ALLOCATABLE :: iobsp 145 !! 146 !! 147 148 ALLOCATE( iobsp(kno) ) 149 150 iobsp(:)=kobsp(:) 151 152 WHERE( iobsp(:) == -1 ) 153 iobsp(:) = 9999999 154 END WHERE 155 156 iobsp(:)=-1*iobsp(:) 157 158 CALL obs_mpp_max_integer( iobsp, kno ) 159 160 kobsp(:)=-1*iobsp(:) 161 162 isum=0 157 163 DO ji = 1, kno 158 IF ( ( kobsi(ji) >= 1 ) .AND. ( kobsi(ji) <= jpi ) .AND. & 159 & ( kobsj(ji) >= 1 ) .AND. ( kobsj(ji) <= jpj ) ) THEN 160 iobsig(ji) = mig( kobsi(ji) ) 161 iobsjg(ji) = mjg( kobsj(ji) ) 162 ELSE 163 iobsig(ji) = -1 164 iobsjg(ji) = -1 164 IF ( kobsp(ji) == 9999999 ) THEN 165 isum=isum+1 166 kobsp(ji)=-1 165 167 ENDIF 166 END DO 167 !----------------------------------------------------------------------- 168 ! Get the decisions from all processors 169 !----------------------------------------------------------------------- 170 ALLOCATE( iobsp(kno,size) ) 171 ALLOCATE( iobsi(kno,size) ) 172 ALLOCATE( iobsj(kno,size) ) 173 CALL mpi_allgather( kobsp, kno, mpi_integer, & 174 & iobsp, kno, mpi_integer, & 175 & mpi_comm_opa, ierr ) 176 CALL mpi_allgather( iobsig, kno, mpi_integer, & 177 & iobsi, kno, mpi_integer, & 178 & mpi_comm_opa, ierr ) 179 CALL mpi_allgather( iobsjg, kno, mpi_integer, & 180 & iobsj, kno, mpi_integer, & 181 & mpi_comm_opa, ierr ) 182 183 !----------------------------------------------------------------------- 184 ! Find the processor with observations from the lowest processor 185 ! number among processors holding the observation. 186 !----------------------------------------------------------------------- 187 kobsp(:) = -1 188 num_sus_obs = 0 189 DO ji = 1, kno 190 DO jj = 1, size 191 IF ( ( kobsp(ji) == -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 192 kobsp(ji) = iobsp(ji,jj) 193 iobsip = iobsi(ji,jj) 194 iobsjp = iobsj(ji,jj) 195 ENDIF 196 IF ( ( kobsp(ji) /= -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 197 IF ( ( iobsip /= iobsi(ji,jj) ) .OR. & 198 & ( iobsjp /= iobsj(ji,jj) ) ) THEN 199 IF ( ( kobsp(ji) < 1000000 ) .AND. & 200 & ( iobsp(ji,jj) < 1000000 ) ) THEN 201 num_sus_obs=num_sus_obs+1 202 ENDIF 203 ENDIF 204 IF ( mppmap(iobsip,iobsjp) /= ( kobsp(ji)+1 ) ) THEN 205 IF ( ( iobsi(ji,jj) /= -1 ) .AND. & 206 & ( iobsj(ji,jj) /= -1 ) ) THEN 207 IF ((mppmap(iobsi(ji,jj),iobsj(ji,jj)) == (iobsp(ji,jj)+1))& 208 & .OR. ( iobsp(ji,jj) < kobsp(ji) ) ) THEN 209 kobsp(ji) = iobsp(ji,jj) 210 iobsip = iobsi(ji,jj) 211 iobsjp = iobsj(ji,jj) 212 ENDIF 213 ENDIF 214 ENDIF 215 ENDIF 216 END DO 217 END DO 218 IF (lwp) WRITE(numout,*) 'Number of suspicious observations: ',num_sus_obs 219 220 DEALLOCATE( iobsj ) 221 DEALLOCATE( iobsi ) 168 ENDDO 169 170 171 IF ( isum > 0 ) THEN 172 IF (lwp) WRITE(numout,*) isum, ' observations failed the grid search.' 173 IF (lwp) WRITE(numout,*)'If ln_grid_search_lookup=.TRUE., try reducing grid_search_res' 174 ENDIF 175 222 176 DEALLOCATE( iobsp ) 177 223 178 #else 224 179 ! no MPI: empty routine 225 #endif 226 !180 #endif 181 227 182 END SUBROUTINE obs_mpp_find_obs_proc 228 183 -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r6486 r10120 7 7 8 8 !!---------------------------------------------------------------------- 9 !! obs_pro_opt : Compute the model counterpart of temperature and 10 !! salinity observations from profiles 11 !! obs_sla_opt : Compute the model counterpart of sea level anomaly 12 !! observations 13 !! obs_sst_opt : Compute the model counterpart of sea surface temperature 14 !! observations 15 !! obs_sss_opt : Compute the model counterpart of sea surface salinity 16 !! observations 17 !! obs_seaice_opt : Compute the model counterpart of sea ice concentration 18 !! observations 19 !! 20 !! obs_vel_opt : Compute the model counterpart of zonal and meridional 21 !! components of velocity from observations. 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 !! obs_surf_opt : Compute the model counterpart of surface data 22 11 !!---------------------------------------------------------------------- 23 12 24 !! * Modules used 13 !! * Modules used 25 14 USE par_kind, ONLY : & ! Precision variables 26 15 & wp 27 16 USE in_out_manager ! I/O manager 28 17 USE obs_inter_sup ! Interpolation support 29 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs ervationpt18 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 30 19 & obs_int_h2d, & 31 20 & obs_int_h2d_init 32 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 33 26 & obs_int_z1d, & 34 27 & obs_int_z1d_spl 35 USE obs_const, ONLY : &36 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 37 30 USE dom_oce, ONLY : & 38 & glamt, glam u, glamv, &39 & gphit, gphi u, gphiv40 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 41 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 37 USE obs_grid, ONLY : & 38 & obs_level_search 42 39 43 40 IMPLICIT NONE … … 46 43 PRIVATE 47 44 48 PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations 49 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 50 & obs_sst_opt, & ! Compute the model counterpart of SST observations 51 & obs_sss_opt, & ! Compute the model counterpart of SSS observations 52 & obs_seaice_opt, & 53 & obs_vel_opt ! Compute the model counterpart of velocity profile data 54 55 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_surf_opt ! Compute the model counterpart of surface obs 47 48 INTEGER, PARAMETER, PUBLIC :: & 49 & imaxavtypes = 20 ! Max number of daily avgd obs types 56 50 57 51 !!---------------------------------------------------------------------- … … 61 55 !!---------------------------------------------------------------------- 62 56 57 !! * Substitutions 58 # include "domzgr_substitute.h90" 63 59 CONTAINS 64 60 65 SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 66 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 67 & kdailyavtypes ) 61 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 63 & kit000, kdaystp, & 64 & pvar1, pvar2, pgdept, pgdepw, & 65 & pmask1, pmask2, & 66 & plam1, plam2, pphi1, pphi2, & 67 & k1dint, k2dint, kdailyavtypes ) 68 68 69 !!----------------------------------------------------------------------- 69 70 !! … … 78 79 !! 79 80 !! First, a vertical profile of horizontally interpolated model 80 !! now temperatures is computed at the obs (lon, lat) point.81 !! now values is computed at the obs (lon, lat) point. 81 82 !! Several horizontal interpolation schemes are available: 82 83 !! - distance-weighted (great circle) (k2dint = 0) … … 86 87 !! - polynomial (quadrilateral grid) (k2dint = 4) 87 88 !! 88 !! Next, the vertical temperatureprofile is interpolated to the89 !! Next, the vertical profile is interpolated to the 89 90 !! data depth points. Two vertical interpolation schemes are 90 91 !! available: … … 96 97 !! routine. 97 98 !! 98 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is99 !! If the logical is switched on, the model equivalent is 99 100 !! a daily mean model temperature field. So, we first compute 100 101 !! the mean, then interpolate only at the end of the day. 101 102 !! 102 !! Note: thein situ temperature observations must be converted103 !! Note: in situ temperature observations must be converted 103 104 !! to potential temperature (the model variable) prior to 104 105 !! assimilation. 105 !!??????????????????????????????????????????????????????????????106 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???107 !!??????????????????????????????????????????????????????????????108 106 !! 109 107 !! ** Action : … … 115 113 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 116 114 !! ! 07-03 (K. Mogensen) General handling of profiles 115 !! ! 15-02 (M. Martin) Combined routine for all profile types 116 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 117 117 !!----------------------------------------------------------------------- 118 118 119 119 !! * Modules used 120 120 USE obs_profiles_def ! Definition of storage space for profile obs. … … 123 123 124 124 !! * Arguments 125 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 126 INTEGER, INTENT(IN) :: kt ! Time step 127 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 125 TYPE(obs_prof), INTENT(INOUT) :: & 126 & prodatqc ! Subset of profile data passing QC 127 INTEGER, INTENT(IN) :: kt ! Time step 128 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 128 129 INTEGER, INTENT(IN) :: kpj 129 130 INTEGER, INTENT(IN) :: kpk 130 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step131 132 INTEGER, INTENT(IN) :: k1dint 133 INTEGER, INTENT(IN) :: k2dint 134 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day131 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 132 ! (kit000-1 = restart time) 133 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 134 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 135 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 135 136 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 136 & ptn, & ! Model temperature field 137 & psn, & ! Model salinity field 138 & ptmask ! Land-sea mask 139 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 140 & pgdept ! Model array of depth levels 137 & pvar1, & ! Model field 1 138 & pvar2, & ! Model field 2 139 & pmask1, & ! Land-sea mask 1 140 & pmask2 ! Land-sea mask 2 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 142 & plam1, & ! Model longitudes for variable 1 143 & plam2, & ! Model longitudes for variable 2 144 & pphi1, & ! Model latitudes for variable 1 145 & pphi2 ! Model latitudes for variable 2 146 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 147 & pgdept, & ! Model array of depth T levels 148 & pgdepw ! Model array of depth W levels 141 149 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 142 & kdailyavtypes! Types for daily averages 150 & kdailyavtypes ! Types for daily averages 151 143 152 !! * Local declarations 144 153 INTEGER :: ji … … 152 161 INTEGER :: iend 153 162 INTEGER :: iobs 163 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 164 INTEGER :: inum_obs 154 165 INTEGER, DIMENSION(imaxavtypes) :: & 155 166 & idailyavtypes 167 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 168 & igrdi1, & 169 & igrdi2, & 170 & igrdj1, & 171 & igrdj2 172 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 173 156 174 REAL(KIND=wp) :: zlam 157 175 REAL(KIND=wp) :: zphi 158 176 REAL(KIND=wp) :: zdaystp 159 177 REAL(KIND=wp), DIMENSION(kpk) :: & 160 & zobsmask, & 178 & zobsmask1, & 179 & zobsmask2, & 161 180 & zobsk, & 162 181 & zobs2k 163 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 182 REAL(KIND=wp), DIMENSION(2,2,1) :: & 183 & zweig1, & 184 & zweig2, & 164 185 & zweig 165 186 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 166 & zmask, & 167 & zintt, & 168 & zints, & 169 & zinmt, & 170 & zinms 187 & zmask1, & 188 & zmask2, & 189 & zint1, & 190 & zint2, & 191 & zinm1, & 192 & zinm2, & 193 & zgdept, & 194 & zgdepw 171 195 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 172 & zglam, & 173 & zgphi 174 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 175 & igrdi, & 176 & igrdj 196 & zglam1, & 197 & zglam2, & 198 & zgphi1, & 199 & zgphi2 200 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 201 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 202 203 LOGICAL :: ld_dailyav 177 204 178 205 !------------------------------------------------------------------------ 179 206 ! Local initialization 180 207 !------------------------------------------------------------------------ 181 ! ...Record and data counters208 ! Record and data counters 182 209 inrc = kt - kit000 + 2 183 210 ipro = prodatqc%npstp(inrc) 184 211 185 212 ! Daily average types 213 ld_dailyav = .FALSE. 186 214 IF ( PRESENT(kdailyavtypes) ) THEN 187 215 idailyavtypes(:) = kdailyavtypes(:) 216 IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 188 217 ELSE 189 218 idailyavtypes(:) = -1 190 219 ENDIF 191 220 192 ! Initialize daily mean for first timestep 221 ! Daily means are calculated for values over timesteps: 222 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 193 223 idayend = MOD( kt - kit000 + 1, kdaystp ) 194 224 195 ! Added kt == 0 test to catch restart case 196 IF ( idayend == 1 .OR. kt == 0) THEN 197 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 225 IF ( ld_dailyav ) THEN 226 227 ! Initialize daily mean for first timestep of the day 228 IF ( idayend == 1 .OR. kt == 0 ) THEN 229 DO jk = 1, jpk 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 prodatqc%vdmean(ji,jj,jk,1) = 0.0 233 prodatqc%vdmean(ji,jj,jk,2) = 0.0 234 END DO 235 END DO 236 END DO 237 ENDIF 238 198 239 DO jk = 1, jpk 199 240 DO jj = 1, jpj 200 241 DO ji = 1, jpi 201 prodatqc%vdmean(ji,jj,jk,1) = 0.0 202 prodatqc%vdmean(ji,jj,jk,2) = 0.0 242 ! Increment field 1 for computing daily mean 243 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 244 & + pvar1(ji,jj,jk) 245 ! Increment field 2 for computing daily mean 246 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 247 & + pvar2(ji,jj,jk) 203 248 END DO 204 249 END DO 205 250 END DO 206 ENDIF 207 208 DO jk = 1, jpk 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 ! Increment the temperature field for computing daily mean 212 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 213 & + ptn(ji,jj,jk) 214 ! Increment the salinity field for computing daily mean 215 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 216 & + psn(ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 ! Compute the daily mean at the end of day 222 zdaystp = 1.0 / REAL( kdaystp ) 223 IF ( idayend == 0 ) THEN 224 DO jk = 1, jpk 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 228 & * zdaystp 229 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 230 & * zdaystp 251 252 ! Compute the daily mean at the end of day 253 zdaystp = 1.0 / REAL( kdaystp ) 254 IF ( idayend == 0 ) THEN 255 IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 256 CALL FLUSH(numout) 257 DO jk = 1, jpk 258 DO jj = 1, jpj 259 DO ji = 1, jpi 260 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 261 & * zdaystp 262 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 263 & * zdaystp 264 END DO 231 265 END DO 232 266 END DO 233 END DO 267 ENDIF 268 234 269 ENDIF 235 270 236 271 ! Get the data for interpolation 237 272 ALLOCATE( & 238 & igrdi(2,2,ipro), & 239 & igrdj(2,2,ipro), & 240 & zglam(2,2,ipro), & 241 & zgphi(2,2,ipro), & 242 & zmask(2,2,kpk,ipro), & 243 & zintt(2,2,kpk,ipro), & 244 & zints(2,2,kpk,ipro) & 273 & igrdi1(2,2,ipro), & 274 & igrdi2(2,2,ipro), & 275 & igrdj1(2,2,ipro), & 276 & igrdj2(2,2,ipro), & 277 & zglam1(2,2,ipro), & 278 & zglam2(2,2,ipro), & 279 & zgphi1(2,2,ipro), & 280 & zgphi2(2,2,ipro), & 281 & zmask1(2,2,kpk,ipro), & 282 & zmask2(2,2,kpk,ipro), & 283 & zint1(2,2,kpk,ipro), & 284 & zint2(2,2,kpk,ipro), & 285 & zgdept(2,2,kpk,ipro), & 286 & zgdepw(2,2,kpk,ipro) & 245 287 & ) 246 288 247 289 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 248 290 iobs = jobs - prodatqc%nprofup 249 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 250 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 251 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 252 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 253 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 254 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 255 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 256 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 291 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 292 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 293 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 294 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 295 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 296 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 297 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 298 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 299 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 300 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 301 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 302 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 303 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 304 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 305 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 306 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 257 307 END DO 258 308 259 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 260 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 261 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 262 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 263 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 309 ! Initialise depth arrays 310 zgdept(:,:,:,:) = 0.0 311 zgdepw(:,:,:,:) = 0.0 312 313 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 314 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 315 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 316 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 317 318 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 319 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 320 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 322 323 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 324 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 264 325 265 326 ! At the end of the day also get interpolated means 266 IF ( idayend == 0 ) THEN327 IF ( ld_dailyav .AND. idayend == 0 ) THEN 267 328 268 329 ALLOCATE( & 269 & zinm t(2,2,kpk,ipro), &270 & zinm s(2,2,kpk,ipro) &330 & zinm1(2,2,kpk,ipro), & 331 & zinm2(2,2,kpk,ipro) & 271 332 & ) 272 333 273 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &274 & prodatqc%vdmean(:,:,:,1), zinm t)275 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &276 & prodatqc%vdmean(:,:,:,2), zinm s)334 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 335 & prodatqc%vdmean(:,:,:,1), zinm1 ) 336 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 337 & prodatqc%vdmean(:,:,:,2), zinm2 ) 277 338 278 339 ENDIF 279 340 341 ! Return if no observations to process 342 ! Has to be done after comm commands to ensure processors 343 ! stay in sync 344 IF ( ipro == 0 ) RETURN 345 280 346 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 281 347 … … 283 349 284 350 IF ( kt /= prodatqc%mstp(jobs) ) THEN 285 351 286 352 IF(lwp) THEN 287 353 WRITE(numout,*) … … 298 364 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 299 365 ENDIF 300 366 301 367 zlam = prodatqc%rlam(jobs) 302 368 zphi = prodatqc%rphi(jobs) 369 370 ! Horizontal weights 371 ! Masked values are calculated later. 372 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 373 374 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 375 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 376 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 377 378 ENDIF 379 380 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 381 382 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 383 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 384 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 385 386 ENDIF 387 388 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 389 390 zobsk(:) = obfillflt 391 392 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 393 394 IF ( idayend == 0 ) THEN 395 ! Daily averaged data 396 397 ! vertically interpolate all 4 corners 398 ista = prodatqc%npvsta(jobs,1) 399 iend = prodatqc%npvend(jobs,1) 400 inum_obs = iend - ista + 1 401 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 402 403 DO iin=1,2 404 DO ijn=1,2 405 406 IF ( k1dint == 1 ) THEN 407 CALL obs_int_z1d_spl( kpk, & 408 & zinm1(iin,ijn,:,iobs), & 409 & zobs2k, zgdept(iin,ijn,:,iobs), & 410 & zmask1(iin,ijn,:,iobs)) 411 ENDIF 412 413 CALL obs_level_search(kpk, & 414 & zgdept(iin,ijn,:,iobs), & 415 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 416 & iv_indic) 417 418 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 419 & prodatqc%var(1)%vdep(ista:iend), & 420 & zinm1(iin,ijn,:,iobs), & 421 & zobs2k, interp_corner(iin,ijn,:), & 422 & zgdept(iin,ijn,:,iobs), & 423 & zmask1(iin,ijn,:,iobs)) 424 425 ENDDO 426 ENDDO 427 428 ENDIF !idayend 429 430 ELSE 431 432 ! Point data 433 434 ! vertically interpolate all 4 corners 435 ista = prodatqc%npvsta(jobs,1) 436 iend = prodatqc%npvend(jobs,1) 437 inum_obs = iend - ista + 1 438 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 439 DO iin=1,2 440 DO ijn=1,2 441 442 IF ( k1dint == 1 ) THEN 443 CALL obs_int_z1d_spl( kpk, & 444 & zint1(iin,ijn,:,iobs),& 445 & zobs2k, zgdept(iin,ijn,:,iobs), & 446 & zmask1(iin,ijn,:,iobs)) 447 448 ENDIF 449 450 CALL obs_level_search(kpk, & 451 & zgdept(iin,ijn,:,iobs),& 452 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 453 & iv_indic) 454 455 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 456 & prodatqc%var(1)%vdep(ista:iend), & 457 & zint1(iin,ijn,:,iobs), & 458 & zobs2k,interp_corner(iin,ijn,:), & 459 & zgdept(iin,ijn,:,iobs), & 460 & zmask1(iin,ijn,:,iobs) ) 303 461 304 ! Horizontal weights and vertical mask 305 306 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 307 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 308 309 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 310 & zglam(:,:,iobs), zgphi(:,:,iobs), & 311 & zmask(:,:,:,iobs), zweig, zobsmask ) 312 313 ENDIF 314 315 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 462 ENDDO 463 ENDDO 464 465 ENDIF 466 467 !------------------------------------------------------------- 468 ! Compute the horizontal interpolation for every profile level 469 !------------------------------------------------------------- 470 471 DO ikn=1,inum_obs 472 iend=ista+ikn-1 473 474 zweig(:,:,1) = 0._wp 475 476 ! This code forces the horizontal weights to be 477 ! zero IF the observation is below the bottom of the 478 ! corners of the interpolation nodes, Or if it is in 479 ! the mask. This is important for observations near 480 ! steep bathymetry 481 DO iin=1,2 482 DO ijn=1,2 483 484 depth_loop1: DO ik=kpk,2,-1 485 IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN 486 487 zweig(iin,ijn,1) = & 488 & zweig1(iin,ijn,1) * & 489 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 490 & - prodatqc%var(1)%vdep(iend)),0._wp) 491 492 EXIT depth_loop1 493 494 ENDIF 495 496 ENDDO depth_loop1 497 498 ENDDO 499 ENDDO 500 501 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 502 & prodatqc%var(1)%vmod(iend:iend) ) 503 504 ! Set QC flag for any observations found below the bottom 505 ! needed as the check here is more strict than that in obs_prep 506 IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 507 508 ENDDO 509 510 DEALLOCATE(interp_corner,iv_indic) 511 512 ENDIF 513 514 ! For the second variable 515 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 316 516 317 517 zobsk(:) = obfillflt 318 518 319 519 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 320 520 321 521 IF ( idayend == 0 ) THEN 522 ! Daily averaged data 523 524 ! vertically interpolate all 4 corners 525 ista = prodatqc%npvsta(jobs,2) 526 iend = prodatqc%npvend(jobs,2) 527 inum_obs = iend - ista + 1 528 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 529 530 DO iin=1,2 531 DO ijn=1,2 532 533 IF ( k1dint == 1 ) THEN 534 CALL obs_int_z1d_spl( kpk, & 535 & zinm2(iin,ijn,:,iobs), & 536 & zobs2k, zgdept(iin,ijn,:,iobs), & 537 & zmask2(iin,ijn,:,iobs)) 538 ENDIF 539 540 CALL obs_level_search(kpk, & 541 & zgdept(iin,ijn,:,iobs), & 542 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 543 & iv_indic) 544 545 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 546 & prodatqc%var(2)%vdep(ista:iend), & 547 & zinm2(iin,ijn,:,iobs), & 548 & zobs2k, interp_corner(iin,ijn,:), & 549 & zgdept(iin,ijn,:,iobs), & 550 & zmask2(iin,ijn,:,iobs)) 551 552 ENDDO 553 ENDDO 554 555 ENDIF !idayend 556 557 ELSE 558 559 ! Point data 560 561 ! vertically interpolate all 4 corners 562 ista = prodatqc%npvsta(jobs,2) 563 iend = prodatqc%npvend(jobs,2) 564 inum_obs = iend - ista + 1 565 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 566 DO iin=1,2 567 DO ijn=1,2 568 569 IF ( k1dint == 1 ) THEN 570 CALL obs_int_z1d_spl( kpk, & 571 & zint2(iin,ijn,:,iobs),& 572 & zobs2k, zgdept(iin,ijn,:,iobs), & 573 & zmask2(iin,ijn,:,iobs)) 574 575 ENDIF 576 577 CALL obs_level_search(kpk, & 578 & zgdept(iin,ijn,:,iobs),& 579 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 580 & iv_indic) 581 582 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 583 & prodatqc%var(2)%vdep(ista:iend), & 584 & zint2(iin,ijn,:,iobs), & 585 & zobs2k,interp_corner(iin,ijn,:), & 586 & zgdept(iin,ijn,:,iobs), & 587 & zmask2(iin,ijn,:,iobs) ) 588 589 ENDDO 590 ENDDO 591 592 ENDIF 593 594 !------------------------------------------------------------- 595 ! Compute the horizontal interpolation for every profile level 596 !------------------------------------------------------------- 597 598 DO ikn=1,inum_obs 599 iend=ista+ikn-1 322 600 323 ! Daily averaged moored buoy (MRB) data 324 325 CALL obs_int_h2d( kpk, kpk, & 326 & zweig, zinmt(:,:,:,iobs), zobsk ) 327 328 329 ELSE 330 331 CALL ctl_stop( ' A nonzero' // & 332 & ' number of profile T BUOY data should' // & 333 & ' only occur at the end of a given day' ) 334 335 ENDIF 336 337 ELSE 338 339 ! Point data 340 341 CALL obs_int_h2d( kpk, kpk, & 342 & zweig, zintt(:,:,:,iobs), zobsk ) 343 344 ENDIF 345 346 !------------------------------------------------------------- 347 ! Compute vertical second-derivative of the interpolating 348 ! polynomial at obs points 349 !------------------------------------------------------------- 350 351 IF ( k1dint == 1 ) THEN 352 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 353 & pgdept, zobsmask ) 354 ENDIF 355 356 !----------------------------------------------------------------- 357 ! Vertical interpolation to the observation point 358 !----------------------------------------------------------------- 359 ista = prodatqc%npvsta(jobs,1) 360 iend = prodatqc%npvend(jobs,1) 361 CALL obs_int_z1d( kpk, & 362 & prodatqc%var(1)%mvk(ista:iend), & 363 & k1dint, iend - ista + 1, & 364 & prodatqc%var(1)%vdep(ista:iend), & 365 & zobsk, zobs2k, & 366 & prodatqc%var(1)%vmod(ista:iend), & 367 & pgdept, zobsmask ) 368 369 ENDIF 370 371 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 372 373 zobsk(:) = obfillflt 374 375 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 376 377 IF ( idayend == 0 ) THEN 378 379 ! Daily averaged moored buoy (MRB) data 380 381 CALL obs_int_h2d( kpk, kpk, & 382 & zweig, zinms(:,:,:,iobs), zobsk ) 383 384 ELSE 385 386 CALL ctl_stop( ' A nonzero' // & 387 & ' number of profile S BUOY data should' // & 388 & ' only occur at the end of a given day' ) 389 390 ENDIF 391 392 ELSE 393 394 ! Point data 395 396 CALL obs_int_h2d( kpk, kpk, & 397 & zweig, zints(:,:,:,iobs), zobsk ) 398 399 ENDIF 400 401 402 !------------------------------------------------------------- 403 ! Compute vertical second-derivative of the interpolating 404 ! polynomial at obs points 405 !------------------------------------------------------------- 406 407 IF ( k1dint == 1 ) THEN 408 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 409 & pgdept, zobsmask ) 410 ENDIF 411 412 !---------------------------------------------------------------- 413 ! Vertical interpolation to the observation point 414 !---------------------------------------------------------------- 415 ista = prodatqc%npvsta(jobs,2) 416 iend = prodatqc%npvend(jobs,2) 417 CALL obs_int_z1d( kpk, & 418 & prodatqc%var(2)%mvk(ista:iend),& 419 & k1dint, iend - ista + 1, & 420 & prodatqc%var(2)%vdep(ista:iend),& 421 & zobsk, zobs2k, & 422 & prodatqc%var(2)%vmod(ista:iend),& 423 & pgdept, zobsmask ) 424 425 ENDIF 426 427 END DO 601 zweig(:,:,1) = 0._wp 602 603 ! This code forces the horizontal weights to be 604 ! zero IF the observation is below the bottom of the 605 ! corners of the interpolation nodes, Or if it is in 606 ! the mask. This is important for observations near 607 ! steep bathymetry 608 DO iin=1,2 609 DO ijn=1,2 610 611 depth_loop2: DO ik=kpk,2,-1 612 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 613 614 zweig(iin,ijn,1) = & 615 & zweig2(iin,ijn,1) * & 616 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 617 & - prodatqc%var(2)%vdep(iend)),0._wp) 618 619 EXIT depth_loop2 620 621 ENDIF 622 623 ENDDO depth_loop2 624 625 ENDDO 626 ENDDO 627 628 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 629 & prodatqc%var(2)%vmod(iend:iend) ) 630 631 ! Set QC flag for any observations found below the bottom 632 ! needed as the check here is more strict than that in obs_prep 633 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 428 634 635 ENDDO 636 637 DEALLOCATE(interp_corner,iv_indic) 638 639 ENDIF 640 641 ENDDO 642 429 643 ! Deallocate the data for interpolation 430 644 DEALLOCATE( & 431 & igrdi, & 432 & igrdj, & 433 & zglam, & 434 & zgphi, & 435 & zmask, & 436 & zintt, & 437 & zints & 645 & igrdi1, & 646 & igrdi2, & 647 & igrdj1, & 648 & igrdj2, & 649 & zglam1, & 650 & zglam2, & 651 & zgphi1, & 652 & zgphi2, & 653 & zmask1, & 654 & zmask2, & 655 & zint1, & 656 & zint2, & 657 & zgdept, & 658 & zgdepw & 438 659 & ) 660 439 661 ! At the end of the day also get interpolated means 440 IF ( idayend == 0 ) THEN662 IF ( ld_dailyav .AND. idayend == 0 ) THEN 441 663 DEALLOCATE( & 442 & zinm t, &443 & zinm s&664 & zinm1, & 665 & zinm2 & 444 666 & ) 445 667 ENDIF 446 668 447 669 prodatqc%nprofup = prodatqc%nprofup + ipro 448 449 END SUBROUTINE obs_pro_opt 450 451 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 & psshn, psshmask, k2dint ) 670 671 END SUBROUTINE obs_prof_opt 672 673 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 674 & kit000, kdaystp, psurf, psurfmask, & 675 & k2dint, ldnightav, plamscl, pphiscl, & 676 & lindegrees ) 677 453 678 !!----------------------------------------------------------------------- 454 679 !! 455 !! *** ROUTINE obs_s la_opt ***456 !! 457 !! ** Purpose : Compute the model counterpart of s ea level anomaly680 !! *** ROUTINE obs_surf_opt *** 681 !! 682 !! ** Purpose : Compute the model counterpart of surface 458 683 !! data by interpolating from the model grid to the 459 684 !! observation point. … … 462 687 !! the model values at the corners of the surrounding grid box. 463 688 !! 464 !! The n ow model SSHis first computed at the obs (lon, lat) point.689 !! The new model value is first computed at the obs (lon, lat) point. 465 690 !! 466 691 !! Several horizontal interpolation schemes are available: … … 470 695 !! - bilinear (quadrilateral grid) (k2dint = 3) 471 696 !! - polynomial (quadrilateral grid) (k2dint = 4) 472 !! 473 !! The sea level anomaly at the observation points is then computed 474 !! by removing a mean dynamic topography (defined at the obs. point). 697 !! 698 !! Two horizontal averaging schemes are also available: 699 !! - weighted radial footprint (k2dint = 5) 700 !! - weighted rectangular footprint (k2dint = 6) 701 !! 475 702 !! 476 703 !! ** Action : … … 478 705 !! History : 479 706 !! ! 07-03 (A. Weaver) 707 !! ! 15-02 (M. Martin) Combined routine for surface types 708 !! ! 17-03 (M. Martin) Added horizontal averaging options 480 709 !!----------------------------------------------------------------------- 481 710 482 711 !! * Modules used 483 712 USE obs_surf_def ! Definition of storage space for surface observations … … 486 715 487 716 !! * Arguments 488 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening 489 INTEGER, INTENT(IN) :: kt ! Time step 490 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 717 TYPE(obs_surf), INTENT(INOUT) :: & 718 & surfdataqc ! Subset of surface data passing QC 719 INTEGER, INTENT(IN) :: kt ! Time step 720 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 491 721 INTEGER, INTENT(IN) :: kpj 492 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 493 ! (kit000-1 = restart time) 494 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 495 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 496 & psshn, & ! Model SSH field 497 & psshmask ! Land-sea mask 498 722 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 723 ! (kit000-1 = restart time) 724 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 725 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 726 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 727 & psurf, & ! Model surface field 728 & psurfmask ! Land-sea mask 729 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 730 REAL(KIND=wp), INTENT(IN) :: & 731 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 732 & pphiscl ! This is the full width (rather than half-width) 733 LOGICAL, INTENT(IN) :: & 734 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 735 499 736 !! * Local declarations 500 737 INTEGER :: ji … … 502 739 INTEGER :: jobs 503 740 INTEGER :: inrc 504 INTEGER :: is la741 INTEGER :: isurf 505 742 INTEGER :: iobs 506 REAL(KIND=wp) :: zlam 507 REAL(KIND=wp) :: zphi 508 REAL(KIND=wp) :: zext(1), zobsmask(1) 509 REAL(kind=wp), DIMENSION(2,2,1) :: & 510 & zweig 743 INTEGER :: imaxifp, imaxjfp 744 INTEGER :: imodi, imodj 745 INTEGER :: idayend 746 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 747 & igrdi, & 748 & igrdj, & 749 & igrdip1, & 750 & igrdjp1 751 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 752 & icount_night, & 753 & imask_night 754 REAL(wp) :: zlam 755 REAL(wp) :: zphi 756 REAL(wp), DIMENSION(1) :: zext, zobsmask 757 REAL(wp) :: zdaystp 511 758 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 512 & zmask, & 513 & zsshl, & 514 & zglam, & 515 & zgphi 516 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 517 & igrdi, & 518 & igrdj 759 & zweig, & 760 & zmask, & 761 & zsurf, & 762 & zsurfm, & 763 & zsurftmp, & 764 & zglam, & 765 & zgphi, & 766 & zglamf, & 767 & zgphif 768 769 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 770 & zintmp, & 771 & zouttmp, & 772 & zmeanday ! to compute model sst in region of 24h daylight (pole) 519 773 520 774 !------------------------------------------------------------------------ 521 775 ! Local initialization 522 776 !------------------------------------------------------------------------ 523 ! ...Record and data counters777 ! Record and data counters 524 778 inrc = kt - kit000 + 2 525 isla = sladatqc%nsstp(inrc) 779 isurf = surfdataqc%nsstp(inrc) 780 781 ! Work out the maximum footprint size for the 782 ! interpolation/averaging in model grid-points - has to be even. 783 784 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 785 786 787 IF ( ldnightav ) THEN 788 789 ! Initialize array for night mean 790 IF ( kt == 0 ) THEN 791 ALLOCATE ( icount_night(kpi,kpj) ) 792 ALLOCATE ( imask_night(kpi,kpj) ) 793 ALLOCATE ( zintmp(kpi,kpj) ) 794 ALLOCATE ( zouttmp(kpi,kpj) ) 795 ALLOCATE ( zmeanday(kpi,kpj) ) 796 nday_qsr = -1 ! initialisation flag for nbc_dcy 797 ENDIF 798 799 ! Night-time means are calculated for night-time values over timesteps: 800 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 801 idayend = MOD( kt - kit000 + 1, kdaystp ) 802 803 ! Initialize night-time mean for first timestep of the day 804 IF ( idayend == 1 .OR. kt == 0 ) THEN 805 DO jj = 1, jpj 806 DO ji = 1, jpi 807 surfdataqc%vdmean(ji,jj) = 0.0 808 zmeanday(ji,jj) = 0.0 809 icount_night(ji,jj) = 0 810 END DO 811 END DO 812 ENDIF 813 814 zintmp(:,:) = 0.0 815 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 816 imask_night(:,:) = INT( zouttmp(:,:) ) 817 818 DO jj = 1, jpj 819 DO ji = 1, jpi 820 ! Increment the temperature field for computing night mean and counter 821 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 822 & + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 823 zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj) 824 icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 825 END DO 826 END DO 827 828 ! Compute the night-time mean at the end of the day 829 zdaystp = 1.0 / REAL( kdaystp ) 830 IF ( idayend == 0 ) THEN 831 IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 832 DO jj = 1, jpj 833 DO ji = 1, jpi 834 ! Test if "no night" point 835 IF ( icount_night(ji,jj) > 0 ) THEN 836 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 837 & / REAL( icount_night(ji,jj) ) 838 ELSE 839 !At locations where there is no night (e.g. poles), 840 ! calculate daily mean instead of night-time mean. 841 surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 842 ENDIF 843 END DO 844 END DO 845 ENDIF 846 847 ENDIF 526 848 527 849 ! Get the data for interpolation 528 850 529 851 ALLOCATE( & 530 & igrdi(2,2,isla), & 531 & igrdj(2,2,isla), & 532 & zglam(2,2,isla), & 533 & zgphi(2,2,isla), & 534 & zmask(2,2,isla), & 535 & zsshl(2,2,isla) & 852 & zweig(imaxifp,imaxjfp,1), & 853 & igrdi(imaxifp,imaxjfp,isurf), & 854 & igrdj(imaxifp,imaxjfp,isurf), & 855 & zglam(imaxifp,imaxjfp,isurf), & 856 & zgphi(imaxifp,imaxjfp,isurf), & 857 & zmask(imaxifp,imaxjfp,isurf), & 858 & zsurf(imaxifp,imaxjfp,isurf), & 859 & zsurftmp(imaxifp,imaxjfp,isurf), & 860 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 861 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 862 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 863 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 536 864 & ) 537 538 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 539 iobs = jobs - sladatqc%nsurfup 540 igrdi(1,1,iobs) = sladatqc%mi(jobs)-1 541 igrdj(1,1,iobs) = sladatqc%mj(jobs)-1 542 igrdi(1,2,iobs) = sladatqc%mi(jobs)-1 543 igrdj(1,2,iobs) = sladatqc%mj(jobs) 544 igrdi(2,1,iobs) = sladatqc%mi(jobs) 545 igrdj(2,1,iobs) = sladatqc%mj(jobs)-1 546 igrdi(2,2,iobs) = sladatqc%mi(jobs) 547 igrdj(2,2,iobs) = sladatqc%mj(jobs) 865 866 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 867 iobs = jobs - surfdataqc%nsurfup 868 DO ji = 0, imaxifp 869 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 870 871 !Deal with wrap around in longitude 872 IF ( imodi < 1 ) imodi = imodi + jpiglo 873 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 874 875 DO jj = 0, imaxjfp 876 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 877 !If model values are out of the domain to the north/south then 878 !set them to be the edge of the domain 879 IF ( imodj < 1 ) imodj = 1 880 IF ( imodj > jpjglo ) imodj = jpjglo 881 882 igrdip1(ji+1,jj+1,iobs) = imodi 883 igrdjp1(ji+1,jj+1,iobs) = imodj 884 885 IF ( ji >= 1 .AND. jj >= 1 ) THEN 886 igrdi(ji,jj,iobs) = imodi 887 igrdj(ji,jj,iobs) = imodj 888 ENDIF 889 890 END DO 891 END DO 548 892 END DO 549 893 550 CALL obs_int_comm_2d( 2, 2, isla, &894 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 551 895 & igrdi, igrdj, glamt, zglam ) 552 CALL obs_int_comm_2d( 2, 2, isla, &896 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 553 897 & igrdi, igrdj, gphit, zgphi ) 554 CALL obs_int_comm_2d( 2, 2, isla, & 555 & igrdi, igrdj, psshmask, zmask ) 556 CALL obs_int_comm_2d( 2, 2, isla, & 557 & igrdi, igrdj, psshn, zsshl ) 898 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 899 & igrdi, igrdj, psurfmask, zmask ) 900 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 901 & igrdi, igrdj, psurf, zsurf ) 902 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 903 & igrdip1, igrdjp1, glamf, zglamf ) 904 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 905 & igrdip1, igrdjp1, gphif, zgphif ) 906 907 ! At the end of the day get interpolated means 908 IF ( idayend == 0 .AND. ldnightav ) THEN 909 910 ALLOCATE( & 911 & zsurfm(imaxifp,imaxjfp,isurf) & 912 & ) 913 914 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 915 & surfdataqc%vdmean(:,:), zsurfm ) 916 917 ENDIF 558 918 559 919 ! Loop over observations 560 561 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 562 563 iobs = jobs - sladatqc%nsurfup 564 565 IF ( kt /= sladatqc%mstp(jobs) ) THEN 566 920 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 921 922 iobs = jobs - surfdataqc%nsurfup 923 924 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 925 567 926 IF(lwp) THEN 568 927 WRITE(numout,*) … … 574 933 WRITE(numout,*) ' Record = ', jobs, & 575 934 & ' kt = ', kt, & 576 & ' mstp = ', s ladatqc%mstp(jobs), &577 & ' ntyp = ', s ladatqc%ntyp(jobs)935 & ' mstp = ', surfdataqc%mstp(jobs), & 936 & ' ntyp = ', surfdataqc%ntyp(jobs) 578 937 ENDIF 579 CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 580 938 CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 939 940 ENDIF 941 942 zlam = surfdataqc%rlam(jobs) 943 zphi = surfdataqc%rphi(jobs) 944 945 IF ( ldnightav .AND. idayend == 0 ) THEN 946 ! Night-time averaged data 947 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 948 ELSE 949 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 950 ENDIF 951 952 IF ( k2dint <= 4 ) THEN 953 954 ! Get weights to interpolate the model value to the observation point 955 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 956 & zglam(:,:,iobs), zgphi(:,:,iobs), & 957 & zmask(:,:,iobs), zweig, zobsmask ) 958 959 ! Interpolate the model value to the observation point 960 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 961 962 ELSE 963 964 ! Get weights to average the model SLA to the observation footprint 965 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 966 & zglam(:,:,iobs), zgphi(:,:,iobs), & 967 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 968 & zmask(:,:,iobs), plamscl, pphiscl, & 969 & lindegrees, zweig, zobsmask ) 970 971 ! Average the model SST to the observation footprint 972 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 973 & zweig, zsurftmp(:,:,iobs), zext ) 974 975 ENDIF 976 977 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 978 ! ... Remove the MDT from the SSH at the observation point to get the SLA 979 surfdataqc%rext(jobs,1) = zext(1) 980 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 981 ELSE 982 surfdataqc%rmod(jobs,1) = zext(1) 581 983 ENDIF 582 984 583 zlam = sladatqc%rlam(jobs) 584 zphi = sladatqc%rphi(jobs) 585 586 ! Get weights to interpolate the model SSH to the observation point 587 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 588 & zglam(:,:,iobs), zgphi(:,:,iobs), & 589 & zmask(:,:,iobs), zweig, zobsmask ) 590 591 592 ! Interpolate the model SSH to the observation point 593 CALL obs_int_h2d( 1, 1, & 594 & zweig, zsshl(:,:,iobs), zext ) 595 596 sladatqc%rext(jobs,1) = zext(1) 597 ! ... Remove the MDT at the observation point 598 sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 985 IF ( zext(1) == obfillflt ) THEN 986 ! If the observation value is a fill value, set QC flag to bad 987 surfdataqc%nqc(jobs) = 4 988 ENDIF 599 989 600 990 END DO … … 602 992 ! Deallocate the data for interpolation 603 993 DEALLOCATE( & 994 & zweig, & 604 995 & igrdi, & 605 996 & igrdj, & … … 607 998 & zgphi, & 608 999 & zmask, & 609 & zsshl & 1000 & zsurf, & 1001 & zsurftmp, & 1002 & zglamf, & 1003 & zgphif, & 1004 & igrdip1,& 1005 & igrdjp1 & 610 1006 & ) 611 1007 612 sladatqc%nsurfup = sladatqc%nsurfup + isla 613 614 END SUBROUTINE obs_sla_opt 615 616 SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 617 & psstn, psstmask, k2dint, ld_nightav ) 618 !!----------------------------------------------------------------------- 619 !! 620 !! *** ROUTINE obs_sst_opt *** 621 !! 622 !! ** Purpose : Compute the model counterpart of surface temperature 623 !! data by interpolating from the model grid to the 624 !! observation point. 625 !! 626 !! ** Method : Linearly interpolate to each observation point using 627 !! the model values at the corners of the surrounding grid box. 628 !! 629 !! The now model SST is first computed at the obs (lon, lat) point. 630 !! 631 !! Several horizontal interpolation schemes are available: 632 !! - distance-weighted (great circle) (k2dint = 0) 633 !! - distance-weighted (small angle) (k2dint = 1) 634 !! - bilinear (geographical grid) (k2dint = 2) 635 !! - bilinear (quadrilateral grid) (k2dint = 3) 636 !! - polynomial (quadrilateral grid) (k2dint = 4) 637 !! 638 !! 639 !! ** Action : 640 !! 641 !! History : 642 !! ! 07-07 (S. Ricci ) : Original 643 !! 644 !!----------------------------------------------------------------------- 645 646 !! * Modules used 647 USE obs_surf_def ! Definition of storage space for surface observations 648 USE sbcdcy 649 650 IMPLICIT NONE 651 652 !! * Arguments 653 TYPE(obs_surf), INTENT(INOUT) :: & 654 & sstdatqc ! Subset of surface data not failing screening 655 INTEGER, INTENT(IN) :: kt ! Time step 656 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 657 INTEGER, INTENT(IN) :: kpj 658 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 659 ! (kit000-1 = restart time) 660 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 661 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 662 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 663 & psstn, & ! Model SST field 664 & psstmask ! Land-sea mask 665 666 !! * Local declarations 667 INTEGER :: ji 668 INTEGER :: jj 669 INTEGER :: jobs 670 INTEGER :: inrc 671 INTEGER :: isst 672 INTEGER :: iobs 673 INTEGER :: idayend 674 REAL(KIND=wp) :: zlam 675 REAL(KIND=wp) :: zphi 676 REAL(KIND=wp) :: zext(1), zobsmask(1) 677 REAL(KIND=wp) :: zdaystp 678 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 679 & icount_sstnight, & 680 & imask_night 681 REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 682 & zintmp, & 683 & zouttmp, & 684 & zmeanday ! to compute model sst in region of 24h daylight (pole) 685 REAL(kind=wp), DIMENSION(2,2,1) :: & 686 & zweig 687 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 688 & zmask, & 689 & zsstl, & 690 & zsstm, & 691 & zglam, & 692 & zgphi 693 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 694 & igrdi, & 695 & igrdj 696 LOGICAL, INTENT(IN) :: ld_nightav 697 698 !----------------------------------------------------------------------- 699 ! Local initialization 700 !----------------------------------------------------------------------- 701 ! ... Record and data counters 702 inrc = kt - kit000 + 2 703 isst = sstdatqc%nsstp(inrc) 704 705 IF ( ld_nightav ) THEN 706 707 ! Initialize array for night mean 708 709 IF ( kt .EQ. 0 ) THEN 710 ALLOCATE ( icount_sstnight(kpi,kpj) ) 711 ALLOCATE ( imask_night(kpi,kpj) ) 712 ALLOCATE ( zintmp(kpi,kpj) ) 713 ALLOCATE ( zouttmp(kpi,kpj) ) 714 ALLOCATE ( zmeanday(kpi,kpj) ) 715 nday_qsr = -1 ! initialisation flag for nbc_dcy 716 ENDIF 717 718 ! Initialize daily mean for first timestep 719 idayend = MOD( kt - kit000 + 1, kdaystp ) 720 721 ! Added kt == 0 test to catch restart case 722 IF ( idayend == 1 .OR. kt == 0) THEN 723 IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 724 DO jj = 1, jpj 725 DO ji = 1, jpi 726 sstdatqc%vdmean(ji,jj) = 0.0 727 zmeanday(ji,jj) = 0.0 728 icount_sstnight(ji,jj) = 0 729 END DO 730 END DO 731 ENDIF 732 733 zintmp(:,:) = 0.0 734 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 735 imask_night(:,:) = INT( zouttmp(:,:) ) 736 737 DO jj = 1, jpj 738 DO ji = 1, jpi 739 ! Increment the temperature field for computing night mean and counter 740 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 741 & + psstn(ji,jj)*imask_night(ji,jj) 742 zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj) 743 icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 744 END DO 745 END DO 746 747 ! Compute the daily mean at the end of day 748 749 zdaystp = 1.0 / REAL( kdaystp ) 750 751 IF ( idayend == 0 ) THEN 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 ! Test if "no night" point 755 IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 756 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 757 & / icount_sstnight(ji,jj) 758 ELSE 759 sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 760 ENDIF 761 END DO 762 END DO 763 ENDIF 764 765 ENDIF 766 767 ! Get the data for interpolation 768 769 ALLOCATE( & 770 & igrdi(2,2,isst), & 771 & igrdj(2,2,isst), & 772 & zglam(2,2,isst), & 773 & zgphi(2,2,isst), & 774 & zmask(2,2,isst), & 775 & zsstl(2,2,isst) & 776 & ) 777 778 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 779 iobs = jobs - sstdatqc%nsurfup 780 igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 781 igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 782 igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 783 igrdj(1,2,iobs) = sstdatqc%mj(jobs) 784 igrdi(2,1,iobs) = sstdatqc%mi(jobs) 785 igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 786 igrdi(2,2,iobs) = sstdatqc%mi(jobs) 787 igrdj(2,2,iobs) = sstdatqc%mj(jobs) 788 END DO 789 790 CALL obs_int_comm_2d( 2, 2, isst, & 791 & igrdi, igrdj, glamt, zglam ) 792 CALL obs_int_comm_2d( 2, 2, isst, & 793 & igrdi, igrdj, gphit, zgphi ) 794 CALL obs_int_comm_2d( 2, 2, isst, & 795 & igrdi, igrdj, psstmask, zmask ) 796 CALL obs_int_comm_2d( 2, 2, isst, & 797 & igrdi, igrdj, psstn, zsstl ) 798 799 ! At the end of the day get interpolated means 800 IF ( idayend == 0 .AND. ld_nightav ) THEN 801 802 ALLOCATE( & 803 & zsstm(2,2,isst) & 804 & ) 805 806 CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 807 & sstdatqc%vdmean(:,:), zsstm ) 808 809 ENDIF 810 811 ! Loop over observations 812 813 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 814 815 iobs = jobs - sstdatqc%nsurfup 816 817 IF ( kt /= sstdatqc%mstp(jobs) ) THEN 818 819 IF(lwp) THEN 820 WRITE(numout,*) 821 WRITE(numout,*) ' E R R O R : Observation', & 822 & ' time step is not consistent with the', & 823 & ' model time step' 824 WRITE(numout,*) ' =========' 825 WRITE(numout,*) 826 WRITE(numout,*) ' Record = ', jobs, & 827 & ' kt = ', kt, & 828 & ' mstp = ', sstdatqc%mstp(jobs), & 829 & ' ntyp = ', sstdatqc%ntyp(jobs) 830 ENDIF 831 CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 832 833 ENDIF 834 835 zlam = sstdatqc%rlam(jobs) 836 zphi = sstdatqc%rphi(jobs) 837 838 ! Get weights to interpolate the model SST to the observation point 839 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 840 & zglam(:,:,iobs), zgphi(:,:,iobs), & 841 & zmask(:,:,iobs), zweig, zobsmask ) 842 843 ! Interpolate the model SST to the observation point 844 845 IF ( ld_nightav ) THEN 846 847 IF ( idayend == 0 ) THEN 848 ! Daily averaged/diurnal cycle of SST data 849 CALL obs_int_h2d( 1, 1, & 850 & zweig, zsstm(:,:,iobs), zext ) 851 ELSE 852 CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // & 853 & ' number of night SST data should' // & 854 & ' only occur at the end of a given day' ) 855 ENDIF 856 857 ELSE 858 859 CALL obs_int_h2d( 1, 1, & 860 & zweig, zsstl(:,:,iobs), zext ) 861 862 ENDIF 863 sstdatqc%rmod(jobs,1) = zext(1) 864 865 END DO 866 867 ! Deallocate the data for interpolation 868 DEALLOCATE( & 869 & igrdi, & 870 & igrdj, & 871 & zglam, & 872 & zgphi, & 873 & zmask, & 874 & zsstl & 875 & ) 876 877 ! At the end of the day also get interpolated means 878 IF ( idayend == 0 .AND. ld_nightav ) THEN 1008 ! At the end of the day also deallocate night-time mean array 1009 IF ( idayend == 0 .AND. ldnightav ) THEN 879 1010 DEALLOCATE( & 880 & zs stm &1011 & zsurfm & 881 1012 & ) 882 1013 ENDIF 883 884 sstdatqc%nsurfup = sstdatqc%nsurfup + isst 885 886 END SUBROUTINE obs_sst_opt 887 888 SUBROUTINE obs_sss_opt 889 !!----------------------------------------------------------------------- 890 !! 891 !! *** ROUTINE obs_sss_opt *** 892 !! 893 !! ** Purpose : Compute the model counterpart of sea surface salinity 894 !! data by interpolating from the model grid to the 895 !! observation point. 896 !! 897 !! ** Method : 898 !! 899 !! ** Action : 900 !! 901 !! History : 902 !! ! ??-?? 903 !!----------------------------------------------------------------------- 904 905 IMPLICIT NONE 906 907 END SUBROUTINE obs_sss_opt 908 909 SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, & 910 & pseaicen, pseaicemask, k2dint ) 911 912 !!----------------------------------------------------------------------- 913 !! 914 !! *** ROUTINE obs_seaice_opt *** 915 !! 916 !! ** Purpose : Compute the model counterpart of surface temperature 917 !! data by interpolating from the model grid to the 918 !! observation point. 919 !! 920 !! ** Method : Linearly interpolate to each observation point using 921 !! the model values at the corners of the surrounding grid box. 922 !! 923 !! The now model sea ice is first computed at the obs (lon, lat) point. 924 !! 925 !! Several horizontal interpolation schemes are available: 926 !! - distance-weighted (great circle) (k2dint = 0) 927 !! - distance-weighted (small angle) (k2dint = 1) 928 !! - bilinear (geographical grid) (k2dint = 2) 929 !! - bilinear (quadrilateral grid) (k2dint = 3) 930 !! - polynomial (quadrilateral grid) (k2dint = 4) 931 !! 932 !! 933 !! ** Action : 934 !! 935 !! History : 936 !! ! 07-07 (S. Ricci ) : Original 937 !! 938 !!----------------------------------------------------------------------- 939 940 !! * Modules used 941 USE obs_surf_def ! Definition of storage space for surface observations 942 943 IMPLICIT NONE 944 945 !! * Arguments 946 TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening 947 INTEGER, INTENT(IN) :: kt ! Time step 948 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 949 INTEGER, INTENT(IN) :: kpj 950 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 951 ! (kit000-1 = restart time) 952 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 953 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 954 & pseaicen, & ! Model sea ice field 955 & pseaicemask ! Land-sea mask 956 957 !! * Local declarations 958 INTEGER :: ji 959 INTEGER :: jj 960 INTEGER :: jobs 961 INTEGER :: inrc 962 INTEGER :: iseaice 963 INTEGER :: iobs 964 965 REAL(KIND=wp) :: zlam 966 REAL(KIND=wp) :: zphi 967 REAL(KIND=wp) :: zext(1), zobsmask(1) 968 REAL(kind=wp), DIMENSION(2,2,1) :: & 969 & zweig 970 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 971 & zmask, & 972 & zseaicel, & 973 & zglam, & 974 & zgphi 975 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 976 & igrdi, & 977 & igrdj 978 979 !------------------------------------------------------------------------ 980 ! Local initialization 981 !------------------------------------------------------------------------ 982 ! ... Record and data counters 983 inrc = kt - kit000 + 2 984 iseaice = seaicedatqc%nsstp(inrc) 985 986 ! Get the data for interpolation 987 988 ALLOCATE( & 989 & igrdi(2,2,iseaice), & 990 & igrdj(2,2,iseaice), & 991 & zglam(2,2,iseaice), & 992 & zgphi(2,2,iseaice), & 993 & zmask(2,2,iseaice), & 994 & zseaicel(2,2,iseaice) & 995 & ) 996 997 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 998 iobs = jobs - seaicedatqc%nsurfup 999 igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1 1000 igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1 1001 igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1 1002 igrdj(1,2,iobs) = seaicedatqc%mj(jobs) 1003 igrdi(2,1,iobs) = seaicedatqc%mi(jobs) 1004 igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1 1005 igrdi(2,2,iobs) = seaicedatqc%mi(jobs) 1006 igrdj(2,2,iobs) = seaicedatqc%mj(jobs) 1007 END DO 1008 1009 CALL obs_int_comm_2d( 2, 2, iseaice, & 1010 & igrdi, igrdj, glamt, zglam ) 1011 CALL obs_int_comm_2d( 2, 2, iseaice, & 1012 & igrdi, igrdj, gphit, zgphi ) 1013 CALL obs_int_comm_2d( 2, 2, iseaice, & 1014 & igrdi, igrdj, pseaicemask, zmask ) 1015 CALL obs_int_comm_2d( 2, 2, iseaice, & 1016 & igrdi, igrdj, pseaicen, zseaicel ) 1017 1018 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 1019 1020 iobs = jobs - seaicedatqc%nsurfup 1021 1022 IF ( kt /= seaicedatqc%mstp(jobs) ) THEN 1023 1024 IF(lwp) THEN 1025 WRITE(numout,*) 1026 WRITE(numout,*) ' E R R O R : Observation', & 1027 & ' time step is not consistent with the', & 1028 & ' model time step' 1029 WRITE(numout,*) ' =========' 1030 WRITE(numout,*) 1031 WRITE(numout,*) ' Record = ', jobs, & 1032 & ' kt = ', kt, & 1033 & ' mstp = ', seaicedatqc%mstp(jobs), & 1034 & ' ntyp = ', seaicedatqc%ntyp(jobs) 1035 ENDIF 1036 CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' ) 1037 1038 ENDIF 1039 1040 zlam = seaicedatqc%rlam(jobs) 1041 zphi = seaicedatqc%rphi(jobs) 1042 1043 ! Get weights to interpolate the model sea ice to the observation point 1044 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 1045 & zglam(:,:,iobs), zgphi(:,:,iobs), & 1046 & zmask(:,:,iobs), zweig, zobsmask ) 1047 1048 ! ... Interpolate the model sea ice to the observation point 1049 CALL obs_int_h2d( 1, 1, & 1050 & zweig, zseaicel(:,:,iobs), zext ) 1051 1052 seaicedatqc%rmod(jobs,1) = zext(1) 1053 1054 END DO 1055 1056 ! Deallocate the data for interpolation 1057 DEALLOCATE( & 1058 & igrdi, & 1059 & igrdj, & 1060 & zglam, & 1061 & zgphi, & 1062 & zmask, & 1063 & zseaicel & 1064 & ) 1065 1066 seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice 1067 1068 END SUBROUTINE obs_seaice_opt 1069 1070 SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 1071 & pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, & 1072 & ld_dailyav ) 1073 !!----------------------------------------------------------------------- 1074 !! 1075 !! *** ROUTINE obs_vel_opt *** 1076 !! 1077 !! ** Purpose : Compute the model counterpart of velocity profile 1078 !! data by interpolating from the model grid to the 1079 !! observation point. 1080 !! 1081 !! ** Method : Linearly interpolate zonal and meridional components of velocity 1082 !! to each observation point using the model values at the corners of 1083 !! the surrounding grid box. The model velocity components are on a 1084 !! staggered C- grid. 1085 !! 1086 !! For velocity data from the TAO array, the model equivalent is 1087 !! a daily mean velocity field. So, we first compute 1088 !! the mean, then interpolate only at the end of the day. 1089 !! 1090 !! ** Action : 1091 !! 1092 !! History : 1093 !! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles 1094 !! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles 1095 !!----------------------------------------------------------------------- 1096 1097 !! * Modules used 1098 USE obs_profiles_def ! Definition of storage space for profile obs. 1099 1100 IMPLICIT NONE 1101 1102 !! * Arguments 1103 TYPE(obs_prof), INTENT(INOUT) :: & 1104 & prodatqc ! Subset of profile data not failing screening 1105 INTEGER, INTENT(IN) :: kt ! Time step 1106 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1107 INTEGER, INTENT(IN) :: kpj 1108 INTEGER, INTENT(IN) :: kpk 1109 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1110 ! (kit000-1 = restart time) 1111 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 1112 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1113 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 1114 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 1115 & pun, & ! Model zonal component of velocity 1116 & pvn, & ! Model meridional component of velocity 1117 & pumask, & ! Land-sea mask 1118 & pvmask ! Land-sea mask 1119 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 1120 & pgdept ! Model array of depth levels 1121 LOGICAL, INTENT(IN) :: ld_dailyav 1122 1123 !! * Local declarations 1124 INTEGER :: ji 1125 INTEGER :: jj 1126 INTEGER :: jk 1127 INTEGER :: jobs 1128 INTEGER :: inrc 1129 INTEGER :: ipro 1130 INTEGER :: idayend 1131 INTEGER :: ista 1132 INTEGER :: iend 1133 INTEGER :: iobs 1134 INTEGER, DIMENSION(imaxavtypes) :: & 1135 & idailyavtypes 1136 REAL(KIND=wp) :: zlam 1137 REAL(KIND=wp) :: zphi 1138 REAL(KIND=wp) :: zdaystp 1139 REAL(KIND=wp), DIMENSION(kpk) :: & 1140 & zobsmasku, & 1141 & zobsmaskv, & 1142 & zobsmask, & 1143 & zobsk, & 1144 & zobs2k 1145 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 1146 & zweigu,zweigv 1147 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 1148 & zumask, zvmask, & 1149 & zintu, & 1150 & zintv, & 1151 & zinmu, & 1152 & zinmv 1153 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1154 & zglamu, zglamv, & 1155 & zgphiu, zgphiv 1156 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1157 & igrdiu, & 1158 & igrdju, & 1159 & igrdiv, & 1160 & igrdjv 1161 1162 !------------------------------------------------------------------------ 1163 ! Local initialization 1164 !------------------------------------------------------------------------ 1165 ! ... Record and data counters 1166 inrc = kt - kit000 + 2 1167 ipro = prodatqc%npstp(inrc) 1168 1169 ! Initialize daily mean for first timestep 1170 idayend = MOD( kt - kit000 + 1, kdaystp ) 1171 1172 ! Added kt == 0 test to catch restart case 1173 IF ( idayend == 1 .OR. kt == 0) THEN 1174 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 1175 prodatqc%vdmean(:,:,:,1) = 0.0 1176 prodatqc%vdmean(:,:,:,2) = 0.0 1177 ENDIF 1178 1179 ! Increment the zonal velocity field for computing daily mean 1180 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:) 1181 ! Increment the meridional velocity field for computing daily mean 1182 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:) 1183 1184 ! Compute the daily mean at the end of day 1185 zdaystp = 1.0 / REAL( kdaystp ) 1186 IF ( idayend == 0 ) THEN 1187 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp 1188 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp 1189 ENDIF 1190 1191 ! Get the data for interpolation 1192 ALLOCATE( & 1193 & igrdiu(2,2,ipro), & 1194 & igrdju(2,2,ipro), & 1195 & igrdiv(2,2,ipro), & 1196 & igrdjv(2,2,ipro), & 1197 & zglamu(2,2,ipro), zglamv(2,2,ipro), & 1198 & zgphiu(2,2,ipro), zgphiv(2,2,ipro), & 1199 & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), & 1200 & zintu(2,2,kpk,ipro), & 1201 & zintv(2,2,kpk,ipro) & 1202 & ) 1203 1204 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1205 iobs = jobs - prodatqc%nprofup 1206 igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1 1207 igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1 1208 igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1 1209 igrdju(1,2,iobs) = prodatqc%mj(jobs,1) 1210 igrdiu(2,1,iobs) = prodatqc%mi(jobs,1) 1211 igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1 1212 igrdiu(2,2,iobs) = prodatqc%mi(jobs,1) 1213 igrdju(2,2,iobs) = prodatqc%mj(jobs,1) 1214 igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1 1215 igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1 1216 igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1 1217 igrdjv(1,2,iobs) = prodatqc%mj(jobs,2) 1218 igrdiv(2,1,iobs) = prodatqc%mi(jobs,2) 1219 igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1 1220 igrdiv(2,2,iobs) = prodatqc%mi(jobs,2) 1221 igrdjv(2,2,iobs) = prodatqc%mj(jobs,2) 1222 END DO 1223 1224 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu ) 1225 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu ) 1226 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask ) 1227 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu ) 1228 1229 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv ) 1230 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv ) 1231 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask ) 1232 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv ) 1233 1234 ! At the end of the day also get interpolated means 1235 IF ( idayend == 0 ) THEN 1236 1237 ALLOCATE( & 1238 & zinmu(2,2,kpk,ipro), & 1239 & zinmv(2,2,kpk,ipro) & 1240 & ) 1241 1242 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, & 1243 & prodatqc%vdmean(:,:,:,1), zinmu ) 1244 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, & 1245 & prodatqc%vdmean(:,:,:,2), zinmv ) 1246 1247 ENDIF 1248 1249 ! loop over observations 1250 1251 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1252 1253 iobs = jobs - prodatqc%nprofup 1254 1255 IF ( kt /= prodatqc%mstp(jobs) ) THEN 1256 1257 IF(lwp) THEN 1258 WRITE(numout,*) 1259 WRITE(numout,*) ' E R R O R : Observation', & 1260 & ' time step is not consistent with the', & 1261 & ' model time step' 1262 WRITE(numout,*) ' =========' 1263 WRITE(numout,*) 1264 WRITE(numout,*) ' Record = ', jobs, & 1265 & ' kt = ', kt, & 1266 & ' mstp = ', prodatqc%mstp(jobs), & 1267 & ' ntyp = ', prodatqc%ntyp(jobs) 1268 ENDIF 1269 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 1270 ENDIF 1271 1272 zlam = prodatqc%rlam(jobs) 1273 zphi = prodatqc%rphi(jobs) 1274 1275 ! Initialize observation masks 1276 1277 zobsmasku(:) = 0.0 1278 zobsmaskv(:) = 0.0 1279 1280 ! Horizontal weights and vertical mask 1281 1282 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1283 1284 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1285 & zglamu(:,:,iobs), zgphiu(:,:,iobs), & 1286 & zumask(:,:,:,iobs), zweigu, zobsmasku ) 1287 1288 ENDIF 1289 1290 1291 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1292 1293 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1294 & zglamv(:,:,iobs), zgphiv(:,:,iobs), & 1295 & zvmask(:,:,:,iobs), zweigv, zobsmasku ) 1296 1297 ENDIF 1298 1299 ! Ensure that the vertical mask on u and v are consistent. 1300 1301 zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) ) 1302 1303 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1304 1305 zobsk(:) = obfillflt 1306 1307 IF ( ld_dailyav ) THEN 1308 1309 IF ( idayend == 0 ) THEN 1310 1311 ! Daily averaged data 1312 1313 CALL obs_int_h2d( kpk, kpk, & 1314 & zweigu, zinmu(:,:,:,iobs), zobsk ) 1315 1316 1317 ELSE 1318 1319 CALL ctl_stop( ' A nonzero' // & 1320 & ' number of U profile data should' // & 1321 & ' only occur at the end of a given day' ) 1322 1323 ENDIF 1324 1325 ELSE 1326 1327 ! Point data 1328 1329 CALL obs_int_h2d( kpk, kpk, & 1330 & zweigu, zintu(:,:,:,iobs), zobsk ) 1331 1332 ENDIF 1333 1334 !------------------------------------------------------------- 1335 ! Compute vertical second-derivative of the interpolating 1336 ! polynomial at obs points 1337 !------------------------------------------------------------- 1338 1339 IF ( k1dint == 1 ) THEN 1340 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 1341 & pgdept, zobsmask ) 1342 ENDIF 1343 1344 !----------------------------------------------------------------- 1345 ! Vertical interpolation to the observation point 1346 !----------------------------------------------------------------- 1347 ista = prodatqc%npvsta(jobs,1) 1348 iend = prodatqc%npvend(jobs,1) 1349 CALL obs_int_z1d( kpk, & 1350 & prodatqc%var(1)%mvk(ista:iend), & 1351 & k1dint, iend - ista + 1, & 1352 & prodatqc%var(1)%vdep(ista:iend), & 1353 & zobsk, zobs2k, & 1354 & prodatqc%var(1)%vmod(ista:iend), & 1355 & pgdept, zobsmask ) 1356 1357 ENDIF 1358 1359 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1360 1361 zobsk(:) = obfillflt 1362 1363 IF ( ld_dailyav ) THEN 1364 1365 IF ( idayend == 0 ) THEN 1366 1367 ! Daily averaged data 1368 1369 CALL obs_int_h2d( kpk, kpk, & 1370 & zweigv, zinmv(:,:,:,iobs), zobsk ) 1371 1372 ELSE 1373 1374 CALL ctl_stop( ' A nonzero' // & 1375 & ' number of V profile data should' // & 1376 & ' only occur at the end of a given day' ) 1377 1378 ENDIF 1379 1380 ELSE 1381 1382 ! Point data 1383 1384 CALL obs_int_h2d( kpk, kpk, & 1385 & zweigv, zintv(:,:,:,iobs), zobsk ) 1386 1387 ENDIF 1388 1389 1390 !------------------------------------------------------------- 1391 ! Compute vertical second-derivative of the interpolating 1392 ! polynomial at obs points 1393 !------------------------------------------------------------- 1394 1395 IF ( k1dint == 1 ) THEN 1396 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 1397 & pgdept, zobsmask ) 1398 ENDIF 1399 1400 !---------------------------------------------------------------- 1401 ! Vertical interpolation to the observation point 1402 !---------------------------------------------------------------- 1403 ista = prodatqc%npvsta(jobs,2) 1404 iend = prodatqc%npvend(jobs,2) 1405 CALL obs_int_z1d( kpk, & 1406 & prodatqc%var(2)%mvk(ista:iend),& 1407 & k1dint, iend - ista + 1, & 1408 & prodatqc%var(2)%vdep(ista:iend),& 1409 & zobsk, zobs2k, & 1410 & prodatqc%var(2)%vmod(ista:iend),& 1411 & pgdept, zobsmask ) 1412 1413 ENDIF 1414 1415 END DO 1416 1417 ! Deallocate the data for interpolation 1418 DEALLOCATE( & 1419 & igrdiu, & 1420 & igrdju, & 1421 & igrdiv, & 1422 & igrdjv, & 1423 & zglamu, zglamv, & 1424 & zgphiu, zgphiv, & 1425 & zumask, zvmask, & 1426 & zintu, & 1427 & zintv & 1428 & ) 1429 ! At the end of the day also get interpolated means 1430 IF ( idayend == 0 ) THEN 1431 DEALLOCATE( & 1432 & zinmu, & 1433 & zinmv & 1434 & ) 1435 ENDIF 1436 1437 prodatqc%nprofup = prodatqc%nprofup + ipro 1438 1439 END SUBROUTINE obs_vel_opt 1014 1015