Changeset 5659 for branches/2015
- Timestamp:
- 2015-07-31T11:59:15+02:00 (9 years ago)
- Location:
- branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS
- Files:
-
- 1 added
- 19 deleted
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4990 r5659 23 23 USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch 24 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 25 USE obs_read_surf ! Reading and allocation of SLA observations 27 26 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 27 USE obs_prep ! Preparation of obs. (grid search etc). 31 28 USE obs_oper ! Observation operators … … 34 31 USE obs_read_altbias ! Bias treatment for altimeter 35 32 USE obs_profiles_def ! Profile data definitions 36 USE obs_profiles ! Profile data storage37 33 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 34 USE obs_types ! Definitions for observation types 42 35 USE mpp_map ! MPP mapping … … 63 56 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles 64 57 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles 65 LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set66 LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set67 LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles68 58 LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies 69 LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files70 LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files71 59 LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature 72 LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature73 LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data74 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files75 60 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration 76 61 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations 77 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data78 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data79 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data80 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data81 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files82 62 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height 83 63 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity … … 91 71 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS 92 72 73 INTEGER, PUBLIC :: numobtypes !: Number of observation types to read in. 93 74 INTEGER, PUBLIC :: n1dint !: Vertical interpolation method 94 75 INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method 95 76 INTEGER, DIMENSION(:), ALLOCATABLE :: nvarsprof !Number of profile variables 77 INTEGER, DIMENSION(:), ALLOCATABLE :: nextrprof !Number of profile extra variables 78 INTEGER, DIMENSION(:), ALLOCATABLE :: nvarssurf !Number of surface variables 79 INTEGER, DIMENSION(:), ALLOCATABLE :: nextrsurf !Number of surface extra variables 96 80 INTEGER, DIMENSION(imaxavtypes) :: & 97 & endailyavtypes !: ENACT data types which are daily average 98 81 & dailyavtypes !: Data types which are daily average 82 83 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdata ! Initial surface data 84 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdataqc ! Surface data after quality control 85 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdata ! Initial profile data 86 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdataqc ! Profile data after quality control 87 88 CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: obstypesprof 89 CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: obstypessurf 90 91 92 99 93 INTEGER, PARAMETER :: MaxNumFiles = 1000 94 100 95 LOGICAL, DIMENSION(MaxNumFiles) :: & 101 96 & ln_profb_ena, & !: Is the feedback files from ENACT data ? 102 ! !: If so use endailyavtypes97 ! !: If so use dailyavtypes 103 98 & ln_profb_enatim !: Change tim for 820 enact data set. 104 99 … … 106 101 & ln_velfb_av !: Is the velocity feedback files daily average? 107 102 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 108 & ld_enact !: Profile data is ENACT so use endailyavtypes103 & ld_enact !: Profile data is ENACT so use dailyavtypes 109 104 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 110 105 & ld_velav !: Velocity data is daily averaged … … 135 130 !! ! 06-10 (A. Weaver) Cleaning and add controls 136 131 !! ! 07-03 (K. Mogensen) General handling of profiles 132 !! ! 15-02 (M. Martin) Simplification of namelist and code 137 133 !!---------------------------------------------------------------------- 138 134 … … 140 136 141 137 !! * Local declarations 142 CHARACTER(len=128) :: enactfiles(MaxNumFiles)143 CHARACTER(len=128) :: coriofiles(MaxNumFiles)144 138 CHARACTER(len=128) :: profbfiles(MaxNumFiles) 145 CHARACTER(len=128) :: sstfiles(MaxNumFiles) 146 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 147 CHARACTER(len=128) :: slafilesact(MaxNumFiles) 148 CHARACTER(len=128) :: slafilespas(MaxNumFiles) 139 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 149 140 CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 150 CHARACTER(len=128) :: seaicefiles(MaxNumFiles) 151 CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 152 CHARACTER(len=128) :: veladcpfiles(MaxNumFiles) 153 CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 154 CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 155 CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 156 CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 141 CHARACTER(len=128) :: seaicefbfiles(MaxNumFiles) 157 142 CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 158 CHARACTER(LEN=128) :: reysstname159 CHARACTER(LEN=12) :: reysstfmt160 143 CHARACTER(LEN=128) :: bias_file 161 144 CHARACTER(LEN=20) :: datestr=" ", timestr=" " 162 NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, & 163 & ln_sla, ln_sladt, ln_slafb, & 164 & ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, & 165 & enactfiles, coriofiles, profbfiles, & 166 & slafilesact, slafilespas, slafbfiles, & 167 & sstfiles, sstfbfiles, & 168 & ln_seaice, seaicefiles, & 145 146 NAMELIST/namobs/ln_t3d, ln_s3d, ln_sla, ln_sss, ln_ssh, & 147 & ln_sst, ln_seaice, ln_vel3d, & 148 & ln_altbias, ln_nea, ln_grid_global, & 149 & ln_grid_search_lookup, ln_cl4, & 150 & ln_ignmis, ln_s_at_t, ln_sstnight, & 151 & ln_profb_ena, ln_profb_enatim, & 152 & profbfiles, slafbfiles, sssfbfiles, & 153 & sshfbfiles, sstfbfiles, seaicefbfiles, & 154 & velfbfiles, bias_file, grid_search_file, & 169 155 & dobsini, dobsend, n1dint, n2dint, & 170 156 & nmsshc, mdtcorr, mdtcutoff, & 171 & ln_reysst, ln_ghrsst, reysstname, reysstfmt, & 172 & ln_sstnight, & 173 & ln_grid_search_lookup, & 174 & grid_search_file, grid_search_res, & 175 & ln_grid_global, bias_file, ln_altbias, & 176 & endailyavtypes, ln_s_at_t, ln_profb_ena, & 177 & ln_vel3d, ln_velavcur, velavcurfiles, & 178 & ln_velhrcur, velhrcurfiles, & 179 & ln_velavadcp, velavadcpfiles, & 180 & ln_velhradcp, velhradcpfiles, & 181 & ln_velfb, velfbfiles, ln_velfb_av, & 182 & ln_profb_enatim, ln_ignmis, ln_cl4 183 184 INTEGER :: jprofset 185 INTEGER :: jveloset 186 INTEGER :: jvar 187 INTEGER :: jnumenact 188 INTEGER :: jnumcorio 189 INTEGER :: jnumprofb 190 INTEGER :: jnumslaact 191 INTEGER :: jnumslapas 192 INTEGER :: jnumslafb 193 INTEGER :: jnumsst 194 INTEGER :: jnumsstfb 195 INTEGER :: jnumseaice 196 INTEGER :: jnumvelavcur 197 INTEGER :: jnumvelhrcur 198 INTEGER :: jnumvelavadcp 199 INTEGER :: jnumvelhradcp 200 INTEGER :: jnumvelfb 201 INTEGER :: ji 202 INTEGER :: jset 157 & grid_search_res, dailyavtypes 158 159 INTEGER :: jtype 203 160 INTEGER :: ios ! Local integer output status for namelist read 204 LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 205 161 INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilesprof 162 INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilessurf 163 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilesprof 164 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilessurf 165 LOGICAL :: lmask(MaxNumFiles) 206 166 !----------------------------------------------------------------------- 207 167 ! Read namelist parameters 208 168 !----------------------------------------------------------------------- 209 169 210 enactfiles(:) = ''211 coriofiles(:) = ''212 170 profbfiles(:) = '' 213 slafilesact(:) = ''214 slafilespas(:) = ''215 171 slafbfiles(:) = '' 216 sstfiles(:) = ''217 172 sstfbfiles(:) = '' 218 seaicefiles(:) = '' 219 velcurfiles(:) = '' 220 veladcpfiles(:) = '' 221 velavcurfiles(:) = '' 222 velhrcurfiles(:) = '' 223 velavadcpfiles(:) = '' 224 velhradcpfiles(:) = '' 173 seaicefbfiles(:) = '' 225 174 velfbfiles(:) = '' 226 velcurfiles(:) = '' 227 veladcpfiles(:) = '' 228 endailyavtypes(:) = -1 229 endailyavtypes(1) = 820 175 dailyavtypes(:) = -1 176 dailyavtypes(1) = 820 230 177 ln_profb_ena(:) = .FALSE. 231 178 ln_profb_enatim(:) = .TRUE. 232 179 ln_velfb_av(:) = .FALSE. 233 180 ln_ignmis = .FALSE. 234 181 235 182 CALL ini_date( dobsini ) 236 183 CALL fin_date( dobsend ) 237 184 238 185 ! Read Namelist namobs : control observation diagnostics 239 186 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation … … 246 193 IF(lwm) WRITE ( numond, namobs ) 247 194 248 ! Count number of files for each type 249 IF (ln_ena) THEN 250 lmask(:) = .FALSE. 251 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 252 jnumenact = COUNT(lmask) 253 ENDIF 254 IF (ln_cor) THEN 255 lmask(:) = .FALSE. 256 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 257 jnumcorio = COUNT(lmask) 258 ENDIF 259 IF (ln_profb) THEN 260 lmask(:) = .FALSE. 261 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 262 jnumprofb = COUNT(lmask) 263 ENDIF 264 IF (ln_sladt) THEN 265 lmask(:) = .FALSE. 266 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 267 jnumslaact = COUNT(lmask) 268 lmask(:) = .FALSE. 269 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 270 jnumslapas = COUNT(lmask) 271 ENDIF 272 IF (ln_slafb) THEN 273 lmask(:) = .FALSE. 274 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 275 jnumslafb = COUNT(lmask) 276 lmask(:) = .FALSE. 277 ENDIF 278 IF (ln_ghrsst) THEN 279 lmask(:) = .FALSE. 280 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 281 jnumsst = COUNT(lmask) 282 ENDIF 283 IF (ln_sstfb) THEN 284 lmask(:) = .FALSE. 285 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 286 jnumsstfb = COUNT(lmask) 287 lmask(:) = .FALSE. 288 ENDIF 289 IF (ln_seaice) THEN 290 lmask(:) = .FALSE. 291 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 292 jnumseaice = COUNT(lmask) 293 ENDIF 294 IF (ln_velavcur) THEN 295 lmask(:) = .FALSE. 296 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 297 jnumvelavcur = COUNT(lmask) 298 ENDIF 299 IF (ln_velhrcur) THEN 300 lmask(:) = .FALSE. 301 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 302 jnumvelhrcur = COUNT(lmask) 303 ENDIF 304 IF (ln_velavadcp) THEN 305 lmask(:) = .FALSE. 306 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 307 jnumvelavadcp = COUNT(lmask) 308 ENDIF 309 IF (ln_velhradcp) THEN 310 lmask(:) = .FALSE. 311 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 312 jnumvelhradcp = COUNT(lmask) 313 ENDIF 314 IF (ln_velfb) THEN 315 lmask(:) = .FALSE. 316 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 317 jnumvelfb = COUNT(lmask) 318 lmask(:) = .FALSE. 319 ENDIF 320 321 ! Control print 195 !Set up list of observation types to be used 196 numproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 197 numsurftypes = COUNT( (/ln_sla, ln_sss, ln_sst, ln_seaice /) ) 198 IF ( numproftypes > 0 ) THEN 199 200 ALLOCATE( obstypesprof(numproftypes) ) 201 ALLOCATE( jnumfilesprof(numproftypes) ) 202 ALLOCATE( obsfilesprof(numproftypes, MaxNumFiles) ) 203 204 DO jtype = 1, numproftypes 205 IF (ln_t3d .OR. ln_s3d) THEN 206 obsfilesprof(:,jtype) = profbfiles(:) 207 obstypesprof(jtype) = 'prof' 208 ENDIF 209 IF (ln_vel3d) THEN 210 obsfilesprof(:,jtype) = velfbfiles(:) 211 obstypesprof(jtype) = 'vel' 212 ENDIF 213 214 lmask(:) = .FALSE. 215 WHERE (obsfilesprof(jtype,:) /= '') lmask(:) = .TRUE. 216 jnumfilesprof(jtype) = COUNT(lmask) 217 END DO 218 219 ENDIF 220 221 IF ( numsurftypes > 0 ) THEN 222 223 ALLOCATE( obstypessurf(numsurftypes) ) 224 ALLOCATE( jnumfilessurf(numproftypes) ) 225 ALLOCATE( obsfilessurf(numsurftypes, MaxNumFiles) ) 226 227 DO jtype = 1, numsurftypes 228 IF (ln_sla) THEN 229 obsfilessurf(:,jtype) = slafbfiles(:) 230 obstypessurf(jtype) = 'sla' 231 ENDIF 232 IF (ln_sss) THEN 233 obsfilessurf(:,jtype) = sssfbfiles(:) 234 obstypessurf(jtype) = 'sss' 235 ENDIF 236 IF (ln_sst) THEN 237 obsfilessurf(:,jtype) = sstfbfiles(:) 238 obstypessurf(jtype) = 'sst' 239 ENDIF 240 #if defined key_lim2 || defined key_lim3 241 IF (ln_seaice) THEN 242 obsfilessurf(:,jtype) = seaicefbfiles(:) 243 obstypessurf(jtype) = 'seaice' 244 ENDIF 245 #endif 246 247 lmask(:) = .FALSE. 248 WHERE (obsfilessurf(jtype,:) /= '') lmask(:) = .TRUE. 249 jnumfilessurf(jtype) = COUNT(lmask) 250 251 END DO 252 253 ENDIF 254 255 !Write namelist settings to stdout 322 256 IF(lwp) THEN 323 257 WRITE(numout,*) … … 327 261 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 328 262 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 329 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena330 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor331 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb332 263 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 333 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt334 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb335 264 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 336 265 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 337 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst338 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst339 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb340 266 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 341 267 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 342 268 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 343 269 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 344 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 345 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 346 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 347 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 348 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 349 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 270 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global 350 271 WRITE(numout,*) & 351 272 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 352 273 IF (ln_grid_search_lookup) & 353 274 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 354 IF (ln_ena) THEN 355 DO ji = 1, jnumenact 356 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 357 TRIM(enactfiles(ji)) 358 END DO 359 ENDIF 360 IF (ln_cor) THEN 361 DO ji = 1, jnumcorio 362 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 363 TRIM(coriofiles(ji)) 364 END DO 365 ENDIF 366 IF (ln_profb) THEN 367 DO ji = 1, jnumprofb 368 IF (ln_profb_ena(ji)) THEN 369 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 370 TRIM(profbfiles(ji)) 371 ELSE 372 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 373 TRIM(profbfiles(ji)) 374 ENDIF 375 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 376 END DO 377 ENDIF 378 IF (ln_sladt) THEN 379 DO ji = 1, jnumslaact 380 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 381 TRIM(slafilesact(ji)) 382 END DO 383 DO ji = 1, jnumslapas 384 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 385 TRIM(slafilespas(ji)) 386 END DO 387 ENDIF 388 IF (ln_slafb) THEN 389 DO ji = 1, jnumslafb 390 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 391 TRIM(slafbfiles(ji)) 392 END DO 393 ENDIF 394 IF (ln_ghrsst) THEN 395 DO ji = 1, jnumsst 396 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 397 TRIM(sstfiles(ji)) 398 END DO 399 ENDIF 400 IF (ln_sstfb) THEN 401 DO ji = 1, jnumsstfb 402 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 403 TRIM(sstfbfiles(ji)) 404 END DO 405 ENDIF 406 IF (ln_seaice) THEN 407 DO ji = 1, jnumseaice 408 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 409 TRIM(seaicefiles(ji)) 410 END DO 411 ENDIF 412 IF (ln_velavcur) THEN 413 DO ji = 1, jnumvelavcur 414 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 415 TRIM(velavcurfiles(ji)) 416 END DO 417 ENDIF 418 IF (ln_velhrcur) THEN 419 DO ji = 1, jnumvelhrcur 420 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 421 TRIM(velhrcurfiles(ji)) 422 END DO 423 ENDIF 424 IF (ln_velavadcp) THEN 425 DO ji = 1, jnumvelavadcp 426 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 427 TRIM(velavadcpfiles(ji)) 428 END DO 429 ENDIF 430 IF (ln_velhradcp) THEN 431 DO ji = 1, jnumvelhradcp 432 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 433 TRIM(velhradcpfiles(ji)) 434 END DO 435 ENDIF 436 IF (ln_velfb) THEN 437 DO ji = 1, jnumvelfb 438 IF (ln_velfb_av(ji)) THEN 439 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 440 TRIM(velfbfiles(ji)) 441 ELSE 442 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 443 TRIM(velfbfiles(ji)) 444 ENDIF 445 END DO 446 ENDIF 447 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 275 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsin 448 276 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 449 277 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint … … 455 283 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 456 284 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 457 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 285 WRITE(numout,*) ' Daily average types = ', dailyavtypes 286 287 IF ( numproftypes > 0 ) THEN 288 DO jtype = 1, numproftypes 289 DO ji = 1, jnumfilesprof(jtype) 290 WRITE(numout,'(1X,2A)') ' '//obstypesprof(jtype)//' input observation file names = ', & 291 TRIM(obsfilesprof(jtype,ji)) 292 IF ( TRIM(obstypesprof(jtype)) == 'prof' ) & 293 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 294 END DO 295 END DO 296 ENDIF 297 298 IF ( numsurftypes > 0 ) THEN 299 DO jtype = 1, numsurftypes 300 DO ji = 1, jnumfilessurf(jtype) 301 WRITE(numout,'(1X,2A)') ' '//obstypessurf(jtype)//' input observation file names = ', & 302 TRIM(obsfilessurf(jtype,ji)) 303 END DO 304 END DO 305 ENDIF 458 306 459 307 ENDIF … … 470 318 ! Parameter control 471 319 #if defined key_diaobs 472 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 473 & ( .NOT. ln_vel3d ).AND. & 474 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 475 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 320 IF ( numobtypes == 0 ) THEN 476 321 IF(lwp) WRITE(numout,cform_war) 477 322 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & … … 494 339 ! Depending on switches read the various observation types 495 340 !----------------------------------------------------------------------- 496 ! - Temperature/salinity profiles 497 498 IF ( ln_t3d .OR. ln_s3d ) THEN 499 500 ! Set the number of variables for profiles to 2 (T and S) 501 nprofvars = 2 502 ! Set the number of extra variables for profiles to 1 (insitu temp). 503 nprofextr = 1 504 505 ! Count how may insitu data sets we have and allocate data. 506 jprofset = 0 507 IF ( ln_ena ) jprofset = jprofset + 1 508 IF ( ln_cor ) jprofset = jprofset + 1 509 IF ( ln_profb ) jprofset = jprofset + jnumprofb 510 nprofsets = jprofset 511 IF ( nprofsets > 0 ) THEN 512 ALLOCATE(ld_enact(nprofsets)) 513 ALLOCATE(profdata(nprofsets)) 514 ALLOCATE(prodatqc(nprofsets)) 515 ENDIF 516 517 jprofset = 0 518 519 ! ENACT insitu data 520 521 IF ( ln_ena ) THEN 522 523 jprofset = jprofset + 1 341 342 IF ( numproftypes > 0 ) THEN 343 344 ALLOCATE(profdata(numproftypes)) 345 ALLOCATE(profdataqc(numproftypes)) 346 ALLOCATE(nvarsprof(numproftypes)) 347 ALLOCATE(nextrprof(numproftypes)) 524 348 525 ld_enact(jprofset) = .TRUE. 526 527 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 528 & jnumenact, enactfiles(1:jnumenact), & 529 & nprofvars, nprofextr, & 530 & nitend-nit000+2, & 531 & dobsini, dobsend, ln_t3d, ln_s3d, & 532 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 533 & kdailyavtypes = endailyavtypes ) 534 535 DO jvar = 1, 2 536 537 CALL obs_prof_staend( profdata(jprofset), jvar ) 538 349 DO jtype = 1, numproftypes 350 351 nvarsprof(jtype) = 2 352 IF ( TRIM(obstypesprof(jtype)) == 'prof' ) nextrprof(jtype) = 1 353 IF ( TRIM(obstypesprof(jtype)) == 'vel' ) nextrprof(jtype) = 2 354 355 !Read in profile or velocity obs types 356 CALL obs_rea_prof( profdata(jtype), & 357 & jnumfilesprof(jtype), & 358 & obsfilesprof(jtype,1:jnumfilesprof(jtype)), & 359 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 360 & dobsini, dobsend, ln_t3d, ln_s3d, & 361 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 362 & kdailyavtypes = dailyavtypes ) 363 364 DO jvar = 1, nvars 365 CALL obs_prof_staend( profdata(jtype), jvar ) 539 366 END DO 540 541 CALL obs_pre_pro ( profdata(jprofset), prodatqc(jprofset), &367 368 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 542 369 & ln_t3d, ln_s3d, ln_nea, & 543 & kdailyavtypes=endailyavtypes ) 544 545 ENDIF 546 547 ! Coriolis insitu data 548 549 IF ( ln_cor ) THEN 550 551 jprofset = jprofset + 1 552 553 ld_enact(jprofset) = .FALSE. 554 555 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 556 & jnumcorio, coriofiles(1:jnumcorio), & 557 & nprofvars, nprofextr, & 558 & nitend-nit000+2, & 559 & dobsini, dobsend, ln_t3d, ln_s3d, & 560 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 561 562 DO jvar = 1, 2 563 564 CALL obs_prof_staend( profdata(jprofset), jvar ) 565 566 END DO 567 568 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 569 & ln_t3d, ln_s3d, ln_nea ) 570 571 ENDIF 572 573 ! Feedback insitu data 574 575 IF ( ln_profb ) THEN 576 577 DO jset = 1, jnumprofb 578 579 jprofset = jprofset + 1 580 ld_enact (jprofset) = ln_profb_ena(jset) 581 582 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 583 & 1, profbfiles(jset:jset), & 584 & nprofvars, nprofextr, & 585 & nitend-nit000+2, & 586 & dobsini, dobsend, ln_t3d, ln_s3d, & 587 & ln_ignmis, ln_s_at_t, & 588 & ld_enact(jprofset).AND.& 589 & ln_profb_enatim(jset), & 590 & .FALSE., kdailyavtypes = endailyavtypes ) 591 592 DO jvar = 1, 2 593 594 CALL obs_prof_staend( profdata(jprofset), jvar ) 595 596 END DO 597 598 IF ( ld_enact(jprofset) ) THEN 599 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 600 & ln_t3d, ln_s3d, ln_nea, & 601 & kdailyavtypes = endailyavtypes ) 602 ELSE 603 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 604 & ln_t3d, ln_s3d, ln_nea ) 605 ENDIF 606 607 END DO 608 609 ENDIF 610 611 ENDIF 612 613 ! - Sea level anomalies 614 IF ( ln_sla ) THEN 615 ! Set the number of variables for sla to 1 616 nslavars = 1 617 618 ! Set the number of extra variables for sla to 2 619 nslaextr = 2 620 621 ! Set the number of sla data sets to 2 622 nslasets = 0 623 IF ( ln_sladt ) THEN 624 nslasets = nslasets + 2 625 ENDIF 626 IF ( ln_slafb ) THEN 627 nslasets = nslasets + jnumslafb 628 ENDIF 629 630 ALLOCATE(sladata(nslasets)) 631 ALLOCATE(sladatqc(nslasets)) 632 sladata(:)%nsurf=0 633 sladatqc(:)%nsurf=0 634 635 nslasets = 0 636 637 ! AVISO SLA data 638 639 IF ( ln_sladt ) THEN 640 641 ! Active SLA observations 642 643 nslasets = nslasets + 1 644 645 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 646 & slafilesact(1:jnumslaact), & 647 & nslavars, nslaextr, nitend-nit000+2, & 648 & dobsini, dobsend, ln_ignmis, .FALSE. ) 649 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 650 & ln_sla, ln_nea ) 651 652 ! Passive SLA observations 653 654 nslasets = nslasets + 1 655 656 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 657 & slafilespas(1:jnumslapas), & 658 & nslavars, nslaextr, nitend-nit000+2, & 659 & dobsini, dobsend, ln_ignmis, .FALSE. ) 660 661 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 662 & ln_sla, ln_nea ) 663 664 ENDIF 665 666 ! Feedback SLA data 667 668 IF ( ln_slafb ) THEN 669 670 DO jset = 1, jnumslafb 671 672 nslasets = nslasets + 1 673 674 CALL obs_rea_sla( 0, sladata(nslasets), 1, & 675 & slafbfiles(jset:jset), & 676 & nslavars, nslaextr, nitend-nit000+2, & 677 & dobsini, dobsend, ln_ignmis, .FALSE. ) 678 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 679 & ln_sla, ln_nea ) 680 681 END DO 682 683 ENDIF 684 685 CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 686 687 ! read in altimeter bias 688 689 IF ( ln_altbias ) THEN 690 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 691 ENDIF 692 693 ENDIF 694 695 ! - Sea surface height 696 IF ( ln_ssh ) THEN 697 IF(lwp) WRITE(numout,*) ' SSH currently not available' 698 ENDIF 699 700 ! - Sea surface temperature 701 IF ( ln_sst ) THEN 702 703 ! Set the number of variables for sst to 1 704 nsstvars = 1 705 706 ! Set the number of extra variables for sst to 0 707 nsstextr = 0 708 709 nsstsets = 0 710 711 IF (ln_reysst) nsstsets = nsstsets + 1 712 IF (ln_ghrsst) nsstsets = nsstsets + 1 713 IF ( ln_sstfb ) THEN 714 nsstsets = nsstsets + jnumsstfb 715 ENDIF 716 717 ALLOCATE(sstdata(nsstsets)) 718 ALLOCATE(sstdatqc(nsstsets)) 719 ALLOCATE(ld_sstnight(nsstsets)) 720 sstdata(:)%nsurf=0 721 sstdatqc(:)%nsurf=0 722 ld_sstnight(:)=.false. 723 724 nsstsets = 0 725 726 IF (ln_reysst) THEN 727 728 nsstsets = nsstsets + 1 729 730 ld_sstnight(nsstsets) = ln_sstnight 731 732 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 733 & nsstvars, nsstextr, & 734 & nitend-nit000+2, dobsini, dobsend ) 735 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 736 & ln_nea ) 737 738 ENDIF 739 740 IF (ln_ghrsst) THEN 741 742 nsstsets = nsstsets + 1 743 744 ld_sstnight(nsstsets) = ln_sstnight 745 746 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 747 & sstfiles(1:jnumsst), & 748 & nsstvars, nsstextr, nitend-nit000+2, & 749 & dobsini, dobsend, ln_ignmis, .FALSE. ) 750 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 751 & ln_nea ) 752 753 ENDIF 754 755 ! Feedback SST data 756 757 IF ( ln_sstfb ) THEN 758 759 DO jset = 1, jnumsstfb 760 761 nsstsets = nsstsets + 1 762 763 ld_sstnight(nsstsets) = ln_sstnight 764 765 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 766 & sstfbfiles(jset:jset), & 767 & nsstvars, nsstextr, nitend-nit000+2, & 768 & dobsini, dobsend, ln_ignmis, .FALSE. ) 769 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 770 & ln_sst, ln_nea ) 771 772 END DO 773 774 ENDIF 775 776 ENDIF 777 778 ! - Sea surface salinity 779 IF ( ln_sss ) THEN 780 IF(lwp) WRITE(numout,*) ' SSS currently not available' 781 ENDIF 782 783 ! - Sea Ice Concentration 784 785 IF ( ln_seaice ) THEN 786 787 ! Set the number of variables for seaice to 1 788 nseaicevars = 1 789 790 ! Set the number of extra variables for seaice to 0 791 nseaiceextr = 0 792 793 ! Set the number of data sets to 1 794 nseaicesets = 1 795 796 ALLOCATE(seaicedata(nseaicesets)) 797 ALLOCATE(seaicedatqc(nseaicesets)) 798 seaicedata(:)%nsurf=0 799 seaicedatqc(:)%nsurf=0 800 801 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 802 & seaicefiles(1:jnumseaice), & 803 & nseaicevars, nseaiceextr, nitend-nit000+2, & 804 & dobsini, dobsend, ln_ignmis, .FALSE. ) 805 806 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 807 & ln_seaice, ln_nea ) 808 809 ENDIF 810 811 IF (ln_vel3d) THEN 812 813 ! Set the number of variables for profiles to 2 (U and V) 814 nvelovars = 2 815 816 ! Set the number of extra variables for profiles to 2 to store 817 ! rotation parameters 818 nveloextr = 2 819 820 jveloset = 0 821 822 IF ( ln_velavcur ) jveloset = jveloset + 1 823 IF ( ln_velhrcur ) jveloset = jveloset + 1 824 IF ( ln_velavadcp ) jveloset = jveloset + 1 825 IF ( ln_velhradcp ) jveloset = jveloset + 1 826 IF (ln_velfb) jveloset = jveloset + jnumvelfb 827 828 nvelosets = jveloset 829 IF ( nvelosets > 0 ) THEN 830 ALLOCATE( velodata(nvelosets) ) 831 ALLOCATE( veldatqc(nvelosets) ) 832 ALLOCATE( ld_velav(nvelosets) ) 833 ENDIF 834 835 jveloset = 0 836 837 ! Daily averaged data 838 839 IF ( ln_velavcur ) THEN 840 841 jveloset = jveloset + 1 842 843 ld_velav(jveloset) = .TRUE. 844 845 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 846 & velavcurfiles(1:jnumvelavcur), & 847 & nvelovars, nveloextr, & 848 & nitend-nit000+2, & 849 & dobsini, dobsend, ln_ignmis, & 850 & ld_velav(jveloset), & 851 & .FALSE. ) 852 853 DO jvar = 1, 2 854 CALL obs_prof_staend( velodata(jveloset), jvar ) 855 END DO 856 857 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 858 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 859 860 ENDIF 861 862 ! High frequency data 863 864 IF ( ln_velhrcur ) THEN 865 866 jveloset = jveloset + 1 867 868 ld_velav(jveloset) = .FALSE. 869 870 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 871 & velhrcurfiles(1:jnumvelhrcur), & 872 & nvelovars, nveloextr, & 873 & nitend-nit000+2, & 874 & dobsini, dobsend, ln_ignmis, & 875 & ld_velav(jveloset), & 876 & .FALSE. ) 877 878 DO jvar = 1, 2 879 CALL obs_prof_staend( velodata(jveloset), jvar ) 880 END DO 881 882 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 883 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 884 885 ENDIF 886 887 ! Daily averaged data 888 889 IF ( ln_velavadcp ) THEN 890 891 jveloset = jveloset + 1 892 893 ld_velav(jveloset) = .TRUE. 894 895 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 896 & velavadcpfiles(1:jnumvelavadcp), & 897 & nvelovars, nveloextr, & 898 & nitend-nit000+2, & 899 & dobsini, dobsend, ln_ignmis, & 900 & ld_velav(jveloset), & 901 & .FALSE. ) 902 903 DO jvar = 1, 2 904 CALL obs_prof_staend( velodata(jveloset), jvar ) 905 END DO 906 907 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 908 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 909 910 ENDIF 911 912 ! High frequency data 913 914 IF ( ln_velhradcp ) THEN 915 916 jveloset = jveloset + 1 917 918 ld_velav(jveloset) = .FALSE. 919 920 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 921 & velhradcpfiles(1:jnumvelhradcp), & 922 & nvelovars, nveloextr, & 923 & nitend-nit000+2, & 924 & dobsini, dobsend, ln_ignmis, & 925 & ld_velav(jveloset), & 926 & .FALSE. ) 927 928 DO jvar = 1, 2 929 CALL obs_prof_staend( velodata(jveloset), jvar ) 930 END DO 931 932 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 933 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 934 935 ENDIF 936 937 IF ( ln_velfb ) THEN 938 939 DO jset = 1, jnumvelfb 940 941 jveloset = jveloset + 1 942 943 ld_velav(jveloset) = ln_velfb_av(jset) 944 945 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 946 & velfbfiles(jset:jset), & 947 & nvelovars, nveloextr, & 948 & nitend-nit000+2, & 949 & dobsini, dobsend, ln_ignmis, & 950 & ld_velav(jveloset), & 951 & .FALSE. ) 952 953 DO jvar = 1, 2 954 CALL obs_prof_staend( velodata(jveloset), jvar ) 955 END DO 956 957 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 958 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 959 960 961 END DO 962 963 ENDIF 964 965 ENDIF 966 370 & kdailyavtypes = dailyavtypes ) 371 372 END DO 373 374 DEALLOCATE( jnumfilesprof, obsfilesprof ) 375 376 ENDIF 377 378 IF ( numsurftypes > 0 ) THEN 379 380 ALLOCATE(surfdata(numsurftypes)) 381 ALLOCATE(surfdatatqc(numsurftypes)) 382 ALLOCATE(nvarssurf(numsurftypes)) 383 ALLOCATE(nextrsurf(numsurftypes)) 384 385 DO jtype = 1, numsurftypes 386 387 nvarssurf(jtype) = 1 388 nextrsurf(jtype) = 0 389 IF ( TRIM(obstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 390 391 !Read in surface obs types 392 CALL obs_rea_surf( surfdata(jtype), jnumfilessurf(jtype), & 393 & obsfilessurf(jtype,1:jnumfilessurf(jtype)), & 394 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 395 & dobsini, dobsend, ln_ignmis, .FALSE. ) 396 397 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 398 399 IF ( TRIM(obstypessurf(jtype)) == 'sla' ) THEN 400 CALL obs_rea_mdt( surfdataqc(jtype), n2dint ) 401 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), n2dint, bias_file ) 402 ENDIF 403 404 DEALLOCATE( jnumfilessurf, obsfilessurf ) 405 406 END DO 407 967 408 END SUBROUTINE dia_obs_init 968 409 … … 1017 458 !! * Local declarations 1018 459 INTEGER :: idaystp ! Number of timesteps per day 1019 INTEGER :: jprofset ! Profile data set loop variable 1020 INTEGER :: jslaset ! SLA data set loop variable 1021 INTEGER :: jsstset ! SST data set loop variable 1022 INTEGER :: jseaiceset ! sea ice data set loop variable 1023 INTEGER :: jveloset ! velocity profile data loop variable 460 INTEGER :: jtype ! data loop variable 1024 461 INTEGER :: jvar ! Variable number 1025 462 #if ! defined key_lim2 && ! defined key_lim3 … … 1050 487 !----------------------------------------------------------------------- 1051 488 1052 ! - Temperature/salinity profiles 1053 IF ( ln_t3d .OR. ln_s3d ) THEN 1054 DO jprofset = 1, nprofsets 1055 IF ( ld_enact(jprofset) ) THEN 1056 CALL obs_pro_opt( prodatqc(jprofset), & 1057 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1058 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1059 & gdept_1d, tmask, n1dint, n2dint, & 1060 & kdailyavtypes = endailyavtypes ) 1061 ELSE 1062 CALL obs_pro_opt( prodatqc(jprofset), & 1063 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1064 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1065 & gdept_1d, tmask, n1dint, n2dint ) 1066 ENDIF 489 IF ( numproftypes > 0 ) THEN 490 DO jtype = 1, numproftypes 491 492 SELECT CASE ( TRIM(obstypesprof(jtype)) ) 493 CASE('prof') 494 profvar1(:,:,:) = tsn(:,:,:,jp_tem) 495 profvar2(:,:,:) = tsn(:,:,:,jp_sal) 496 profmask1(:,:,:) = tmask(:,:,:) 497 profmask2(:,:,:) = tmask(:,:,:) 498 CASE('vel') 499 profvar1(:,:,:) = un(:,:,:) 500 profvar2(:,:,:) = vn(:,:,:) 501 profmask1(:,:,:) = umask(:,:,:) 502 profmask2(:,:,:) = vmask(:,:,:) 503 END SELECT 504 505 CALL obs_prof_opt( profdataqc(jtype), & 506 & kstp, jpi, jpj, jpk, nit000, idaystp, & 507 & profvar1, profvar2, & 508 & gdept_1d, profmask1, profmask2, n1dint, n2dint, & 509 & kdailyavtypes = dailyavtypes ) 510 1067 511 END DO 1068 ENDIF 1069 1070 ! - Sea surface anomaly 1071 IF ( ln_sla ) THEN 1072 DO jslaset = 1, nslasets 1073 CALL obs_sla_opt( sladatqc(jslaset), & 1074 & kstp, jpi, jpj, nit000, sshn, & 1075 & tmask(:,:,1), n2dint ) 1076 END DO 1077 ENDIF 1078 1079 ! - Sea surface temperature 1080 IF ( ln_sst ) THEN 1081 DO jsstset = 1, nsstsets 1082 CALL obs_sst_opt( sstdatqc(jsstset), & 1083 & kstp, jpi, jpj, nit000, idaystp, & 1084 & tsn(:,:,1,jp_tem), tmask(:,:,1), & 1085 & n2dint, ld_sstnight(jsstset) ) 512 513 ENDIF 514 515 IF ( numsurftypes > 0 ) THEN 516 DO jtype = 1, numsurftypes 517 518 SELECT CASE ( TRIM(obstypessurf(jtype)) ) 519 CASE('sst') 520 surfvar(:,:) = tsn(:,:,1,jp_tem) 521 CASE('sla') 522 surfvar(:,:) = sshn(:,:) 523 CASE('sss') 524 surfvar(:,:) = tsn(:,:,1,jp_sal) 525 #if defined key_lim2 || defined key_lim3 526 CASE('seaice') 527 surfvar(:,:) = 1._wp - frld(:,:) 528 #endif 529 END SELECT 530 531 CALL obs_surf_opt( surfdatqc(jtype), & 532 & kstp, jpi, jpj, nit000, surfvar, & 533 & tmask(:,:,1), n2dint, ld_sstnight ) 534 1086 535 END DO 1087 ENDIF 1088 1089 ! - Sea surface salinity 1090 IF ( ln_sss ) THEN 1091 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1092 ENDIF 1093 1094 #if defined key_lim2 || defined key_lim3 1095 IF ( ln_seaice ) THEN 1096 DO jseaiceset = 1, nseaicesets 1097 CALL obs_seaice_opt( seaicedatqc(jseaiceset), & 1098 & kstp, jpi, jpj, nit000, 1.-frld, & 1099 & tmask(:,:,1), n2dint ) 1100 END DO 1101 ENDIF 1102 #endif 1103 1104 ! - Velocity profiles 1105 IF ( ln_vel3d ) THEN 1106 DO jveloset = 1, nvelosets 1107 ! zonal component of velocity 1108 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, & 1109 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, & 1110 n1dint, n2dint, ld_velav(jveloset) ) 1111 END DO 1112 ENDIF 1113 536 537 ENDIF 538 1114 539 #if ! defined key_lim2 && ! defined key_lim3 1115 540 CALL wrk_dealloc(jpi,jpj,frld) … … 1139 564 !! * Local declarations 1140 565 1141 INTEGER :: jprofset ! Profile data set loop variable 1142 INTEGER :: jveloset ! Velocity data set loop variable 1143 INTEGER :: jslaset ! SLA data set loop variable 1144 INTEGER :: jsstset ! SST data set loop variable 1145 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1146 INTEGER :: jset 1147 INTEGER :: jfbini 1148 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1149 CHARACTER(LEN=10) :: cdtmp 566 INTEGER :: jtype ! Data set loop variable 1150 567 !----------------------------------------------------------------------- 1151 568 ! Depending on switches call various observation output routines 1152 569 !----------------------------------------------------------------------- 1153 570 1154 ! - Temperature/salinity profiles 1155 1156 IF( ln_t3d .OR. ln_s3d ) THEN 1157 1158 ! Copy data from prodatqc to profdata structures 1159 DO jprofset = 1, nprofsets 1160 1161 CALL obs_prof_decompress( prodatqc(jprofset), & 1162 & profdata(jprofset), .TRUE., numout ) 1163 571 IF ( numproftypes > 0 ) THEN 572 DO jtype = 1, numproftypes 573 574 CALL obs_prof_decompress( profdataqc(jtype), & 575 & profdata(jtype), .TRUE., numout ) 576 577 CALL obs_wri_prof( obstypesprof(jtype), profdata(jtype), n2dint ) 578 1164 579 END DO 1165 1166 ! Write the profiles. 1167 1168 jprofset = 0 1169 1170 ! ENACT insitu data 1171 1172 IF ( ln_ena ) THEN 1173 1174 jprofset = jprofset + 1 1175 1176 CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 1177 1178 ENDIF 1179 1180 ! Coriolis insitu data 1181 1182 IF ( ln_cor ) THEN 1183 1184 jprofset = jprofset + 1 1185 1186 CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 1187 1188 ENDIF 1189 1190 ! Feedback insitu data 1191 1192 IF ( ln_profb ) THEN 1193 1194 jfbini = jprofset + 1 1195 1196 DO jprofset = jfbini, nprofsets 1197 1198 jset = jprofset - jfbini + 1 1199 WRITE(cdtmp,'(A,I2.2)')'profb_',jset 1200 CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 1201 1202 END DO 1203 1204 ENDIF 1205 1206 ENDIF 1207 1208 ! - Sea surface anomaly 1209 IF ( ln_sla ) THEN 1210 1211 ! Copy data from sladatqc to sladata structures 1212 DO jslaset = 1, nslasets 1213 1214 CALL obs_surf_decompress( sladatqc(jslaset), & 1215 & sladata(jslaset), .TRUE., numout ) 580 581 ENDIF 582 583 IF ( numsurftypes > 0 ) THEN 584 DO jtype = 1, numsurftypes 585 586 CALL obs_surf_decompress( surfdatqc(jtype), & 587 & surfdata(jtype), .TRUE., numout ) 588 589 CALL obs_wri_surf( obstypessurf(jtype), surfdata(jtype), n2dint ) 1216 590 1217 591 END DO 1218 1219 jslaset = 0 1220 1221 ! Write the AVISO SLA data 1222 1223 IF ( ln_sladt ) THEN 1224 1225 jslaset = 1 1226 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) ) 1227 jslaset = 2 1228 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) ) 1229 1230 ENDIF 1231 1232 IF ( ln_slafb ) THEN 1233 1234 jfbini = jslaset + 1 1235 1236 DO jslaset = jfbini, nslasets 1237 1238 jset = jslaset - jfbini + 1 1239 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset 1240 CALL obs_wri_sla( cdtmp, sladata(jslaset) ) 1241 1242 END DO 1243 1244 ENDIF 1245 1246 ENDIF 1247 1248 ! - Sea surface temperature 1249 IF ( ln_sst ) THEN 1250 1251 ! Copy data from sstdatqc to sstdata structures 1252 DO jsstset = 1, nsstsets 1253 1254 CALL obs_surf_decompress( sstdatqc(jsstset), & 1255 & sstdata(jsstset), .TRUE., numout ) 1256 1257 END DO 1258 1259 jsstset = 0 1260 1261 ! Write the AVISO SST data 1262 1263 IF ( ln_reysst ) THEN 1264 1265 jsstset = jsstset + 1 1266 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) ) 1267 1268 ENDIF 1269 1270 IF ( ln_ghrsst ) THEN 1271 1272 jsstset = jsstset + 1 1273 CALL obs_wri_sst( 'ghr', sstdata(jsstset) ) 1274 1275 ENDIF 1276 1277 IF ( ln_sstfb ) THEN 1278 1279 jfbini = jsstset + 1 1280 1281 DO jsstset = jfbini, nsstsets 1282 1283 jset = jsstset - jfbini + 1 1284 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset 1285 CALL obs_wri_sst( cdtmp, sstdata(jsstset) ) 1286 1287 END DO 1288 1289 ENDIF 1290 1291 ENDIF 1292 1293 ! - Sea surface salinity 1294 IF ( ln_sss ) THEN 1295 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1296 ENDIF 1297 1298 ! - Sea Ice Concentration 1299 IF ( ln_seaice ) THEN 1300 1301 ! Copy data from seaicedatqc to seaicedata structures 1302 DO jseaiceset = 1, nseaicesets 1303 1304 CALL obs_surf_decompress( seaicedatqc(jseaiceset), & 1305 & seaicedata(jseaiceset), .TRUE., numout ) 1306 1307 END DO 1308 1309 ! Write the Sea Ice data 1310 DO jseaiceset = 1, nseaicesets 1311 1312 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset 1313 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) ) 1314 1315 END DO 1316 1317 ENDIF 1318 1319 ! Velocity data 1320 IF( ln_vel3d ) THEN 1321 1322 ! Copy data from veldatqc to velodata structures 1323 DO jveloset = 1, nvelosets 1324 1325 CALL obs_prof_decompress( veldatqc(jveloset), & 1326 & velodata(jveloset), .TRUE., numout ) 1327 1328 END DO 1329 1330 ! Write the profiles. 1331 1332 jveloset = 0 1333 1334 ! Daily averaged data 1335 1336 IF ( ln_velavcur ) THEN 1337 1338 jveloset = jveloset + 1 1339 1340 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint ) 1341 1342 ENDIF 1343 1344 ! High frequency data 1345 1346 IF ( ln_velhrcur ) THEN 1347 1348 jveloset = jveloset + 1 1349 1350 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint ) 1351 1352 ENDIF 1353 1354 ! Daily averaged data 1355 1356 IF ( ln_velavadcp ) THEN 1357 1358 jveloset = jveloset + 1 1359 1360 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint ) 1361 1362 ENDIF 1363 1364 ! High frequency data 1365 1366 IF ( ln_velhradcp ) THEN 1367 1368 jveloset = jveloset + 1 1369 1370 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint ) 1371 1372 ENDIF 1373 1374 ! Feedback velocity data 1375 1376 IF ( ln_velfb ) THEN 1377 1378 jfbini = jveloset + 1 1379 1380 DO jveloset = jfbini, nvelosets 1381 1382 jset = jveloset - jfbini + 1 1383 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset 1384 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint ) 1385 1386 END DO 1387 1388 ENDIF 1389 1390 ENDIF 592 593 ENDIF 594 1391 595 1392 596 END SUBROUTINE dia_obs_wri … … 1409 613 1410 614 !! diaobs deallocation 1411 IF ( nprofsets > 0 ) THEN 1412 DEALLOCATE(ld_enact, & 1413 & profdata, & 1414 & prodatqc) 1415 END IF 1416 IF ( ln_sla ) THEN 1417 DEALLOCATE(sladata, & 1418 & sladatqc) 1419 END IF 1420 IF ( ln_seaice ) THEN 1421 DEALLOCATE(sladata, & 1422 & sladatqc) 1423 END IF 1424 IF ( ln_sst ) THEN 1425 DEALLOCATE(sstdata, & 1426 & sstdatqc) 1427 END IF 1428 IF ( ln_vel3d ) THEN 1429 DEALLOCATE(ld_velav, & 1430 & velodata, & 1431 & veldatqc) 1432 END IF 615 IF ( numproftypes > 0 ) DEALLOCATE(profdata, profdataqc, nvarsprof, nextrprof) 616 IF ( numsurftypes > 0 ) DEALLOCATE(surfdata, surfdataqc, nvarssurf, nextrsurf) 617 1433 618 END SUBROUTINE dia_obs_dealloc 1434 619 -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/julian.F90
r2715 r5659 24 24 & greg2jul ! Convert date to relative time 25 25 26 !! $Id$ 26 27 CONTAINS 27 28 -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r4245 r5659 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 … … 46 35 PRIVATE 47 36 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 37 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile observations 38 & obs_surf_opt ! Compute the model counterpart of surface observations 54 39 55 40 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types … … 63 48 CONTAINS 64 49 65 SUBROUTINE obs_pro _opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &66 & ptn, psn, pgdept, ptmask, k1dint, k2dint, &67 & kdailyavtypes )50 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 51 & pvar1, pvar2, pgdept, pmask1, k1dint, k2dint, & 52 & kdailyavtypes ) 68 53 !!----------------------------------------------------------------------- 69 54 !! … … 78 63 !! 79 64 !! First, a vertical profile of horizontally interpolated model 80 !! now temperatures is computed at the obs (lon, lat) point.65 !! now values is computed at the obs (lon, lat) point. 81 66 !! Several horizontal interpolation schemes are available: 82 67 !! - distance-weighted (great circle) (k2dint = 0) … … 86 71 !! - polynomial (quadrilateral grid) (k2dint = 4) 87 72 !! 88 !! Next, the vertical temperatureprofile is interpolated to the73 !! Next, the vertical profile is interpolated to the 89 74 !! data depth points. Two vertical interpolation schemes are 90 75 !! available: … … 96 81 !! routine. 97 82 !! 98 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is83 !! If the logical is switched on, the model equivalent is 99 84 !! a daily mean model temperature field. So, we first compute 100 85 !! the mean, then interpolate only at the end of the day. 101 86 !! 102 !! Note: thein situ temperature observations must be converted87 !! Note: in situ temperature observations must be converted 103 88 !! to potential temperature (the model variable) prior to 104 89 !! assimilation. 105 !!??????????????????????????????????????????????????????????????106 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???107 !!??????????????????????????????????????????????????????????????108 90 !! 109 91 !! ** Action : … … 115 97 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 116 98 !! ! 07-03 (K. Mogensen) General handling of profiles 99 !! ! 15-02 (M. Martin) Combined routine for all profile types 117 100 !!----------------------------------------------------------------------- 118 101 … … 134 117 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 135 118 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 136 & ptn, & ! Model temperature field 137 & psn, & ! Model salinity field 138 & ptmask ! Land-sea mask 119 & pvar1, & ! Model field 1 120 & pvar2, & ! Model field 2 121 & pmask1, & ! Land-sea mask 1 122 & pmask2 ! Land-sea mask 2 139 123 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 140 124 & pgdept ! Model array of depth levels … … 158 142 REAL(KIND=wp) :: zdaystp 159 143 REAL(KIND=wp), DIMENSION(kpk) :: & 144 & zobsmask1, & 145 & zobsmask2, & 160 146 & zobsmask, & 161 147 & zobsk, & 162 148 & zobs2k 163 149 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 164 & zweig 150 & zweig1, & 151 & zweig2 165 152 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 166 & zmask, & 167 & zintt, & 168 & zints, & 169 & zinmt, & 170 & zinms 153 & zmask1, & 154 & zmask2, & 155 & zint1, & 156 & zint2, & 157 & zinm1, & 158 & zinm2 171 159 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 172 & zglam, & 173 & zgphi 160 & zglam1, & 161 & zglam2, & 162 & zgphi1, & 163 & zgphi2 174 164 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 175 & igrdi, & 176 & igrdj 165 & igrdi1, & 166 & igrdi2, & 167 & igrdj1, & 168 & igrdj2 177 169 178 170 !------------------------------------------------------------------------ … … 209 201 DO jj = 1, jpj 210 202 DO ji = 1, jpi 211 ! Increment the temperature fieldfor computing daily mean203 ! Increment field 1 for computing daily mean 212 204 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 213 & + p tn(ji,jj,jk)214 ! Increment the salinity fieldfor computing daily mean205 & + pvar1(ji,jj,jk) 206 ! Increment field 2 for computing daily mean 215 207 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 216 & + p sn(ji,jj,jk)208 & + pvar2(ji,jj,jk) 217 209 END DO 218 210 END DO … … 236 228 ! Get the data for interpolation 237 229 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) & 230 & igrdi1(2,2,ipro), & 231 & igrdi2(2,2,ipro), & 232 & igrdj1(2,2,ipro), & 233 & igrdj2(2,2,ipro), & 234 & zglam1(2,2,ipro), & 235 & zglam2(2,2,ipro), & 236 & zgphi1(2,2,ipro), & 237 & zgphi2(2,2,ipro), & 238 & zmask1(2,2,kpk,ipro), & 239 & zmask2(2,2,kpk,ipro), & 240 & zint1(2,2,kpk,ipro), & 241 & zint2(2,2,kpk,ipro) & 245 242 & ) 246 243 247 244 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 248 245 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) 246 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 247 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 248 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 249 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 250 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 251 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 252 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 253 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 254 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 255 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 256 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 257 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 258 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 259 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 260 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 261 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 257 262 END DO 258 263 259 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 260 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 261 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 262 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 263 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 264 CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, glamt1, zglam1 ) 265 CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, gphit1, zgphi1 ) 266 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 267 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1, zint1 ) 268 269 CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, glamt2, zglam2 ) 270 CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, gphit2, zgphi2 ) 271 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask1, zmask2 ) 272 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2, zint2 ) 264 273 265 274 ! At the end of the day also get interpolated means … … 267 276 268 277 ALLOCATE( & 269 & zinm t(2,2,kpk,ipro), &270 & zinm s(2,2,kpk,ipro) &278 & zinm1(2,2,kpk,ipro), & 279 & zinm2(2,2,kpk,ipro) & 271 280 & ) 272 281 273 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi , igrdj, &274 & prodatqc%vdmean(:,:,:,1), zinm t)275 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi , igrdj, &276 & prodatqc%vdmean(:,:,:,2), zinm s)282 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, & 283 & prodatqc%vdmean(:,:,:,1), zinm1 ) 284 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, & 285 & prodatqc%vdmean(:,:,:,2), zinm2 ) 277 286 278 287 ENDIF … … 304 313 ! Horizontal weights and vertical mask 305 314 306 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 307 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 315 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 308 316 309 317 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 310 & zglam(:,:,iobs), zgphi(:,:,iobs), & 311 & zmask(:,:,:,iobs), zweig, zobsmask ) 312 318 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 319 & zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 320 321 ENDIF 322 323 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 324 325 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 326 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 327 & zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 328 313 329 ENDIF 314 330 … … 321 337 IF ( idayend == 0 ) THEN 322 338 323 ! Daily averaged moored buoy (MRB)data339 ! Daily averaged data 324 340 325 341 CALL obs_int_h2d( kpk, kpk, & 326 & zweig , zinmt(:,:,:,iobs), zobsk )342 & zweig1, zinm1(:,:,:,iobs), zobsk ) 327 343 328 344 … … 330 346 331 347 CALL ctl_stop( ' A nonzero' // & 332 & ' number of profile T BUOYdata should' // &348 & ' number of average profile data should' // & 333 349 & ' only occur at the end of a given day' ) 334 350 … … 340 356 341 357 CALL obs_int_h2d( kpk, kpk, & 342 & zweig , zintt(:,:,:,iobs), zobsk )358 & zweig1, zint1(:,:,:,iobs), zobsk ) 343 359 344 360 ENDIF … … 377 393 IF ( idayend == 0 ) THEN 378 394 379 ! Daily averaged moored buoy (MRB)data395 ! Daily averaged data 380 396 381 397 CALL obs_int_h2d( kpk, kpk, & 382 & zweig , zinms(:,:,:,iobs), zobsk )398 & zweig2, zinm2(:,:,:,iobs), zobsk ) 383 399 384 400 ELSE 385 401 386 402 CALL ctl_stop( ' A nonzero' // & 387 & ' number of profile S BUOYdata should' // &403 & ' number of average profile data should' // & 388 404 & ' only occur at the end of a given day' ) 389 405 … … 395 411 396 412 CALL obs_int_h2d( kpk, kpk, & 397 & zweig , zints(:,:,:,iobs), zobsk )413 & zweig2, zint2(:,:,:,iobs), zobsk ) 398 414 399 415 ENDIF … … 429 445 ! Deallocate the data for interpolation 430 446 DEALLOCATE( & 431 & igrdi, & 432 & igrdj, & 433 & zglam, & 434 & zgphi, & 435 & zmask, & 436 & zintt, & 437 & zints & 447 & igrdi1, & 448 & igrdi2, & 449 & igrdj1, & 450 & igrdj2, & 451 & zglam1, & 452 & zglam2, & 453 & zgphi1, & 454 & zgphi2, & 455 & zmask1, & 456 & zmask2, & 457 & zint1, & 458 & zint2 & 438 459 & ) 439 460 ! At the end of the day also get interpolated means 440 461 IF ( idayend == 0 ) THEN 441 462 DEALLOCATE( & 442 & zinm t, &443 & zinm s&463 & zinm1, & 464 & zinm2 & 444 465 & ) 445 466 ENDIF … … 447 468 prodatqc%nprofup = prodatqc%nprofup + ipro 448 469 449 END SUBROUTINE obs_pro _opt450 451 SUBROUTINE obs_s la_opt( sladatqc, kt, kpi, kpj, kit000, &452 & ps shn, psshmask, k2dint)470 END SUBROUTINE obs_prof_opt 471 472 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, kit000, & 473 & psurf, psurfmask, k2dint, ld_nightav ) 453 474 !!----------------------------------------------------------------------- 454 475 !! 455 !! *** ROUTINE obs_s la_opt ***456 !! 457 !! ** Purpose : Compute the model counterpart of s ea level anomaly476 !! *** ROUTINE obs_surf_opt *** 477 !! 478 !! ** Purpose : Compute the model counterpart of surface 458 479 !! data by interpolating from the model grid to the 459 480 !! observation point. … … 462 483 !! the model values at the corners of the surrounding grid box. 463 484 !! 464 !! The n ow model SSHis first computed at the obs (lon, lat) point.485 !! The new model value is first computed at the obs (lon, lat) point. 465 486 !! 466 487 !! Several horizontal interpolation schemes are available: … … 471 492 !! - polynomial (quadrilateral grid) (k2dint = 4) 472 493 !! 473 !! The sea level anomaly at the observation points is then computed474 !! by removing a mean dynamic topography (defined at the obs. point).475 494 !! 476 495 !! ** Action : … … 478 497 !! History : 479 498 !! ! 07-03 (A. Weaver) 499 !! ! 15-02 (M. Martin) Combined routine for surface types 480 500 !!----------------------------------------------------------------------- 481 501 … … 486 506 487 507 !! * Arguments 488 TYPE(obs_surf), INTENT(INOUT) :: s ladatqc ! Subset of surface data not failing screening508 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 489 509 INTEGER, INTENT(IN) :: kt ! Time step 490 510 INTEGER, INTENT(IN) :: kpi ! Model grid parameters … … 494 514 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 495 515 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 496 & ps shn, & ! Model SSHfield497 & ps shmask! Land-sea mask516 & psurf, & ! Model surface field 517 & psurfmask ! Land-sea mask 498 518 499 519 !! * Local declarations … … 502 522 INTEGER :: jobs 503 523 INTEGER :: inrc 504 INTEGER :: is la524 INTEGER :: isurf 505 525 INTEGER :: iobs 506 526 REAL(KIND=wp) :: zlam … … 511 531 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 512 532 & zmask, & 513 & zs shl, &533 & zsurf, & 514 534 & zglam, & 515 535 & zgphi … … 523 543 ! ... Record and data counters 524 544 inrc = kt - kit000 + 2 525 isla = sladatqc%nsstp(inrc) 545 isurf = surfdataqc%nsstp(inrc) 546 547 IF ( ld_nightav ) THEN 548 549 ! Initialize array for night mean 550 IF ( kt .EQ. 0 ) THEN 551 ALLOCATE ( icount_night(kpi,kpj) ) 552 ALLOCATE ( imask_night(kpi,kpj) ) 553 ALLOCATE ( zintmp(kpi,kpj) ) 554 ALLOCATE ( zouttmp(kpi,kpj) ) 555 ALLOCATE ( zmeanday(kpi,kpj) ) 556 nday_qsr = -1 ! initialisation flag for nbc_dcy 557 ENDIF 558 559 ! Initialize daily mean for first timestep 560 idayend = MOD( kt - kit000 + 1, kdaystp ) 561 562 ! Added kt == 0 test to catch restart case 563 IF ( idayend == 1 .OR. kt == 0) THEN 564 IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 565 DO jj = 1, jpj 566 DO ji = 1, jpi 567 surfdataqc%vdmean(ji,jj) = 0.0 568 zmeanday(ji,jj) = 0.0 569 icount_night(ji,jj) = 0 570 END DO 571 END DO 572 ENDIF 573 574 zintmp(:,:) = 0.0 575 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 576 imask_night(:,:) = INT( zouttmp(:,:) ) 577 578 DO jj = 1, jpj 579 DO ji = 1, jpi 580 ! Increment the temperature field for computing night mean and counter 581 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 582 & + psurf(ji,jj)*imask_night(ji,jj) 583 zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj) 584 icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 585 END DO 586 END DO 587 588 ! Compute the daily mean at the end of day 589 590 zdaystp = 1.0 / REAL( kdaystp ) 591 592 IF ( idayend == 0 ) THEN 593 DO jj = 1, jpj 594 DO ji = 1, jpi 595 ! Test if "no night" point 596 IF ( icount_night(ji,jj) .NE. 0 ) THEN 597 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 598 & / icount_night(ji,jj) 599 ELSE 600 surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 601 ENDIF 602 END DO 603 END DO 604 ENDIF 605 606 ENDIF 526 607 527 608 ! Get the data for interpolation 528 609 529 610 ALLOCATE( & 530 & igrdi(2,2,is la), &531 & igrdj(2,2,is la), &532 & zglam(2,2,is la), &533 & zgphi(2,2,is la), &534 & zmask(2,2,is la), &535 & zs shl(2,2,isla) &611 & igrdi(2,2,isurf), & 612 & igrdj(2,2,isurf), & 613 & zglam(2,2,isurf), & 614 & zgphi(2,2,isurf), & 615 & zmask(2,2,isurf), & 616 & zsurf(2,2,isurf) & 536 617 & ) 537 618 538 DO jobs = s ladatqc%nsurfup + 1, sladatqc%nsurfup + isla539 iobs = jobs - s ladatqc%nsurfup540 igrdi(1,1,iobs) = s ladatqc%mi(jobs)-1541 igrdj(1,1,iobs) = s ladatqc%mj(jobs)-1542 igrdi(1,2,iobs) = s ladatqc%mi(jobs)-1543 igrdj(1,2,iobs) = s ladatqc%mj(jobs)544 igrdi(2,1,iobs) = s ladatqc%mi(jobs)545 igrdj(2,1,iobs) = s ladatqc%mj(jobs)-1546 igrdi(2,2,iobs) = s ladatqc%mi(jobs)547 igrdj(2,2,iobs) = s ladatqc%mj(jobs)619 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 620 iobs = jobs - surfdataqc%nsurfup 621 igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 622 igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 623 igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 624 igrdj(1,2,iobs) = surfdataqc%mj(jobs) 625 igrdi(2,1,iobs) = surfdataqc%mi(jobs) 626 igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 627 igrdi(2,2,iobs) = surfdataqc%mi(jobs) 628 igrdj(2,2,iobs) = surfdataqc%mj(jobs) 548 629 END DO 549 630 550 CALL obs_int_comm_2d( 2, 2, is la, &631 CALL obs_int_comm_2d( 2, 2, isurf, & 551 632 & igrdi, igrdj, glamt, zglam ) 552 CALL obs_int_comm_2d( 2, 2, is la, &633 CALL obs_int_comm_2d( 2, 2, isurf, & 553 634 & 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 ) 635 CALL obs_int_comm_2d( 2, 2, isurf, & 636 & igrdi, igrdj, psurfmask, zmask ) 637 CALL obs_int_comm_2d( 2, 2, isurf, & 638 & igrdi, igrdj, psurf, zsurf ) 639 640 ! At the end of the day get interpolated means 641 IF ( idayend == 0 .AND. ld_nightav ) THEN 642 643 ALLOCATE( & 644 & zsurfm(2,2,isurf) & 645 & ) 646 647 CALL obs_int_comm_2d( 2, 2, isurf, igrdi, igrdj, & 648 & surfdataqc%vdmean(:,:), zsurfm ) 649 650 ENDIF 558 651 559 652 ! 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 653 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 654 655 iobs = jobs - surfdataqc%nsurfup 656 657 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 566 658 567 659 IF(lwp) THEN … … 574 666 WRITE(numout,*) ' Record = ', jobs, & 575 667 & ' kt = ', kt, & 576 & ' mstp = ', s ladatqc%mstp(jobs), &577 & ' ntyp = ', s ladatqc%ntyp(jobs)668 & ' mstp = ', surfdataqc%mstp(jobs), & 669 & ' ntyp = ', surfdataqc%ntyp(jobs) 578 670 ENDIF 579 671 CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) … … 581 673 ENDIF 582 674 583 zlam = s ladatqc%rlam(jobs)584 zphi = s ladatqc%rphi(jobs)585 586 ! Get weights to interpolate the model SSHto the observation point675 zlam = surfdataqc%rlam(jobs) 676 zphi = surfdataqc%rphi(jobs) 677 678 ! Get weights to interpolate the model value to the observation point 587 679 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 588 680 & zglam(:,:,iobs), zgphi(:,:,iobs), & … … 590 682 591 683 592 ! Interpolate the model SSH to the observation point 593 CALL obs_int_h2d( 1, 1, & 594 & zweig, zsshl(:,:,iobs), zext ) 684 ! Interpolate the model field to the observation point 685 IF ( ld_nightav ) THEN 686 687 IF ( idayend == 0 ) THEN 688 ! Daily averaged data 689 CALL obs_int_h2d( 1, 1, & 690 & zweig, zsurfm(:,:,iobs), zext ) 691 ELSE 692 CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // & 693 & ' number of night data should' // & 694 & ' only occur at the end of a given day' ) 695 ENDIF 696 697 ELSE 698 699 CALL obs_int_h2d( 1, 1, & 700 & zweig, zsurf(:,:,iobs), zext ) 701 702 ENDIF 595 703 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) 704 IF ( surfdataqc%nextr == 2 ) THEN 705 ! ... Remove the MDT from the SSH at the observation point to get the SLA 706 surfdataqc%rext(jobs,1) = zext(1) 707 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 708 ELSE 709 surfdataqc%rmod(jobs,1) = zext(1) 710 ENDIF 599 711 600 712 END DO … … 607 719 & zgphi, & 608 720 & zmask, & 609 & zsshl & 610 & ) 611 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 & 721 & zsurf & 875 722 & ) 876 723 … … 878 725 IF ( idayend == 0 .AND. ld_nightav ) THEN 879 726 DEALLOCATE( & 880 & zs stm &727 & zsurfm & 881 728 & ) 882 729 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 730 731 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 732 733 END SUBROUTINE obs_surf_opt 1440 734 1441 735 END MODULE obs_oper -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r4292 r5659 7 7 8 8 !!--------------------------------------------------------------------- 9 !! obs_pre_pro : First level check and screening of T/S profiles 10 !! obs_pre_sla : First level check and screening of SLA observations 11 !! obs_pre_sst : First level check and screening of SLA observations 12 !! obs_pre_seaice : First level check and screening of sea ice observations 13 !! obs_pre_vel : First level check and screening of velocity obs. 14 !! obs_scr : Basic screening of the observations 15 !! obs_coo_tim : Compute number of time steps to the observation time 16 !! obs_sor : Sort the observation arrays 9 !! obs_pre_prof : First level check and screening of profile observations 10 !! obs_pre_surf : First level check and screening of surface observations 11 !! obs_scr : Basic screening of the observations 12 !! obs_coo_tim : Compute number of time steps to the observation time 13 !! obs_sor : Sort the observation arrays 17 14 !!--------------------------------------------------------------------- 18 15 !! * Modules used … … 36 33 37 34 PUBLIC & 38 & obs_pre_pro, & ! First level check and screening of profiles 39 & obs_pre_sla, & ! First level check and screening of SLA data 40 & obs_pre_sst, & ! First level check and screening of SLA data 41 & obs_pre_seaice, & ! First level check and screening of sea ice data 42 & obs_pre_vel, & ! First level check and screening of velocity profiles 43 & calc_month_len ! Calculate the number of days in the months of a year 35 & obs_pre_prof, & ! First level check and screening of profile obs 36 & obs_pre_surf, & ! First level check and screening of surface obs 37 & calc_month_len ! Calculate the number of days in the months of a year 44 38 45 39 !!---------------------------------------------------------------------- … … 51 45 CONTAINS 52 46 53 SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, & 54 & kdailyavtypes ) 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE obs_pre_pro *** 57 !! 58 !! ** Purpose : First level check and screening of T and S profiles 59 !! 60 !! ** Method : First level check and screening of T and S profiles 61 !! 62 !! ** Action : 63 !! 64 !! References : 65 !! 66 !! History : 67 !! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d 68 !! ! 2007-03 (K. Mogensen) General handling of profiles 69 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 70 !!---------------------------------------------------------------------- 71 !! * Modules used 72 USE domstp ! Domain: set the time-step 73 USE par_oce ! Ocean parameters 74 USE dom_oce, ONLY : & ! Geographical information 75 & glamt, & 76 & gphit, & 77 & gdept_1d,& 78 & tmask, & 79 & nproc 80 !! * Arguments 81 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 82 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 83 LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature 84 LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity 85 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 86 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 87 & kdailyavtypes! Types for daily averages 88 !! * Local declarations 89 INTEGER :: iyea0 ! Initial date 90 INTEGER :: imon0 ! - (year, month, day, hour, minute) 91 INTEGER :: iday0 92 INTEGER :: ihou0 93 INTEGER :: imin0 94 INTEGER :: icycle ! Current assimilation cycle 95 ! Counters for observations that 96 INTEGER :: iotdobs ! - outside time domain 97 INTEGER :: iosdtobs ! - outside space domain (temperature) 98 INTEGER :: iosdsobs ! - outside space domain (salinity) 99 INTEGER :: ilantobs ! - within a model land cell (temperature) 100 INTEGER :: ilansobs ! - within a model land cell (salinity) 101 INTEGER :: inlatobs ! - close to land (temperature) 102 INTEGER :: inlasobs ! - close to land (salinity) 103 INTEGER :: igrdobs ! - fail the grid search 104 ! Global counters for observations that 105 INTEGER :: iotdobsmpp ! - outside time domain 106 INTEGER :: iosdtobsmpp ! - outside space domain (temperature) 107 INTEGER :: iosdsobsmpp ! - outside space domain (salinity) 108 INTEGER :: ilantobsmpp ! - within a model land cell (temperature) 109 INTEGER :: ilansobsmpp ! - within a model land cell (salinity) 110 INTEGER :: inlatobsmpp ! - close to land (temperature) 111 INTEGER :: inlasobsmpp ! - close to land (salinity) 112 INTEGER :: igrdobsmpp ! - fail the grid search 113 TYPE(obs_prof_valid) :: llvalid ! Profile selection 114 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 115 & llvvalid ! T,S selection 116 INTEGER :: jvar ! Variable loop variable 117 INTEGER :: jobs ! Obs. loop variable 118 INTEGER :: jstp ! Time loop variable 119 INTEGER :: inrc ! Time index variable 120 121 IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...' 122 123 ! Initial date initialization (year, month, day, hour, minute) 124 iyea0 = ndate0 / 10000 125 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 126 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 127 ihou0 = 0 128 imin0 = 0 129 130 icycle = no ! Assimilation cycle 131 132 ! Diagnotics counters for various failures. 133 134 iotdobs = 0 135 igrdobs = 0 136 iosdtobs = 0 137 iosdsobs = 0 138 ilantobs = 0 139 ilansobs = 0 140 inlatobs = 0 141 inlasobs = 0 142 143 ! ----------------------------------------------------------------------- 144 ! Find time coordinate for profiles 145 ! ----------------------------------------------------------------------- 146 147 IF ( PRESENT(kdailyavtypes) ) THEN 148 CALL obs_coo_tim_prof( icycle, & 149 & iyea0, imon0, iday0, ihou0, imin0, & 150 & profdata%nprof, profdata%nyea, profdata%nmon, & 151 & profdata%nday, profdata%nhou, profdata%nmin, & 152 & profdata%ntyp, profdata%nqc, profdata%mstp, & 153 & iotdobs, kdailyavtypes = kdailyavtypes ) 154 ELSE 155 CALL obs_coo_tim_prof( icycle, & 156 & iyea0, imon0, iday0, ihou0, imin0, & 157 & profdata%nprof, profdata%nyea, profdata%nmon, & 158 & profdata%nday, profdata%nhou, profdata%nmin, & 159 & profdata%ntyp, profdata%nqc, profdata%mstp, & 160 & iotdobs ) 161 ENDIF 162 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 163 164 ! ----------------------------------------------------------------------- 165 ! Check for profiles failing the grid search 166 ! ----------------------------------------------------------------------- 167 168 CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, & 169 & profdata%nqc, igrdobs ) 170 171 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 172 173 ! ----------------------------------------------------------------------- 174 ! Reject all observations for profiles with nqc > 10 175 ! ----------------------------------------------------------------------- 176 177 CALL obs_pro_rej( profdata ) 178 179 ! ----------------------------------------------------------------------- 180 ! Check for land points. This includes points below the model 181 ! bathymetry so this is done for every point in the profile 182 ! ----------------------------------------------------------------------- 183 184 ! Temperature 185 186 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), & 187 & profdata%npvsta(:,1), profdata%npvend(:,1), & 188 & jpi, jpj, & 189 & jpk, & 190 & profdata%mi, profdata%mj, & 191 & profdata%var(1)%mvk, & 192 & profdata%rlam, profdata%rphi, & 193 & profdata%var(1)%vdep, & 194 & glamt, gphit, & 195 & gdept_1d, tmask, & 196 & profdata%nqc, profdata%var(1)%nvqc, & 197 & iosdtobs, ilantobs, & 198 & inlatobs, ld_nea ) 199 200 CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp ) 201 CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp ) 202 CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp ) 203 204 ! Salinity 205 206 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), & 207 & profdata%npvsta(:,2), profdata%npvend(:,2), & 208 & jpi, jpj, & 209 & jpk, & 210 & profdata%mi, profdata%mj, & 211 & profdata%var(2)%mvk, & 212 & profdata%rlam, profdata%rphi, & 213 & profdata%var(2)%vdep, & 214 & glamt, gphit, & 215 & gdept_1d, tmask, & 216 & profdata%nqc, profdata%var(2)%nvqc, & 217 & iosdsobs, ilansobs, & 218 & inlasobs, ld_nea ) 219 220 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 221 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 222 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 223 224 ! ----------------------------------------------------------------------- 225 ! Copy useful data from the profdata data structure to 226 ! the prodatqc data structure 227 ! ----------------------------------------------------------------------- 228 229 ! Allocate the selection arrays 230 231 ALLOCATE( llvalid%luse(profdata%nprof) ) 232 DO jvar = 1,profdata%nvar 233 ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) ) 234 END DO 235 236 ! We want all data which has qc flags <= 10 237 238 llvalid%luse(:) = ( profdata%nqc(:) <= 10 ) 239 DO jvar = 1,profdata%nvar 240 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 241 END DO 242 243 ! The actual copying 244 245 CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, & 246 & lvalid=llvalid, lvvalid=llvvalid ) 247 248 ! Dellocate the selection arrays 249 DEALLOCATE( llvalid%luse ) 250 DO jvar = 1,profdata%nvar 251 DEALLOCATE( llvvalid(jvar)%luse ) 252 END DO 253 254 ! ----------------------------------------------------------------------- 255 ! Print information about what observations are left after qc 256 ! ----------------------------------------------------------------------- 257 258 ! Update the total observation counter array 259 260 IF(lwp) THEN 261 WRITE(numout,*) 262 WRITE(numout,*) 'obs_pre_pro :' 263 WRITE(numout,*) '~~~~~~~~~~~' 264 WRITE(numout,*) 265 WRITE(numout,*) ' Profiles outside time domain = ', & 266 & iotdobsmpp 267 WRITE(numout,*) ' Remaining profiles that failed grid search = ', & 268 & igrdobsmpp 269 WRITE(numout,*) ' Remaining T data outside space domain = ', & 270 & iosdtobsmpp 271 WRITE(numout,*) ' Remaining T data at land points = ', & 272 & ilantobsmpp 273 IF (ld_nea) THEN 274 WRITE(numout,*) ' Remaining T data near land points (removed) = ',& 275 & inlatobsmpp 276 ELSE 277 WRITE(numout,*) ' Remaining T data near land points (kept) = ',& 278 & inlatobsmpp 279 ENDIF 280 WRITE(numout,*) ' T data accepted = ', & 281 & prodatqc%nvprotmpp(1) 282 WRITE(numout,*) ' Remaining S data outside space domain = ', & 283 & iosdsobsmpp 284 WRITE(numout,*) ' Remaining S data at land points = ', & 285 & ilansobsmpp 286 IF (ld_nea) THEN 287 WRITE(numout,*) ' Remaining S data near land points (removed) = ',& 288 & inlasobsmpp 289 ELSE 290 WRITE(numout,*) ' Remaining S data near land points (kept) = ',& 291 & inlasobsmpp 292 ENDIF 293 WRITE(numout,*) ' S data accepted = ', & 294 & prodatqc%nvprotmpp(2) 295 296 WRITE(numout,*) 297 WRITE(numout,*) ' Number of observations per time step :' 298 WRITE(numout,*) 299 WRITE(numout,997) 300 WRITE(numout,998) 301 ENDIF 302 303 DO jobs = 1, prodatqc%nprof 304 inrc = prodatqc%mstp(jobs) + 2 - nit000 305 prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1 306 DO jvar = 1, prodatqc%nvar 307 IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN 308 prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + & 309 & ( prodatqc%npvend(jobs,jvar) - & 310 & prodatqc%npvsta(jobs,jvar) + 1 ) 311 ENDIF 312 END DO 313 END DO 314 315 316 CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, & 317 & nitend - nit000 + 2 ) 318 DO jvar = 1, prodatqc%nvar 319 CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), & 320 & prodatqc%nvstpmpp(:,jvar), & 321 & nitend - nit000 + 2 ) 322 END DO 323 324 IF ( lwp ) THEN 325 DO jstp = nit000 - 1, nitend 326 inrc = jstp - nit000 + 2 327 WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 328 & prodatqc%nvstpmpp(inrc,1), & 329 & prodatqc%nvstpmpp(inrc,2) 330 END DO 331 ENDIF 332 333 997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity') 334 998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------') 335 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 336 337 END SUBROUTINE obs_pre_pro 338 339 SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea ) 47 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 340 48 !!---------------------------------------------------------------------- 341 49 !! *** ROUTINE obs_pre_sla *** 342 50 !! 343 !! ** Purpose : First level check and screening of SLAobservations344 !! 345 !! ** Method : First level check and screening of SLAobservations51 !! ** Purpose : First level check and screening of surface observations 52 !! 53 !! ** Method : First level check and screening of surface observations 346 54 !! 347 55 !! ** Action : … … 352 60 !! ! 2007-03 (A. Weaver, K. Mogensen) Original 353 61 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 62 !! ! 2015-02 (M. Martin) Combined routine for surface types. 354 63 !!---------------------------------------------------------------------- 355 64 !! * Modules used … … 362 71 & nproc 363 72 !! * Arguments 364 TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data 365 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening 366 LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data 73 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of SLA data 74 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of SLA data not failing screening 367 75 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 368 76 !! * Local declarations … … 411 119 412 120 ! ----------------------------------------------------------------------- 413 ! Find time coordinate for SLAdata121 ! Find time coordinate for surface data 414 122 ! ----------------------------------------------------------------------- 415 123 416 124 CALL obs_coo_tim( icycle, & 417 125 & iyea0, imon0, iday0, ihou0, imin0, & 418 & s ladata%nsurf, sladata%nyea, sladata%nmon, &419 & s ladata%nday, sladata%nhou, sladata%nmin, &420 & s ladata%nqc, sladata%mstp, iotdobs )126 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 127 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 128 & surfdata%nqc, surfdata%mstp, iotdobs ) 421 129 422 130 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 423 131 424 132 ! ----------------------------------------------------------------------- 425 ! Check for SLAdata failing the grid search426 ! ----------------------------------------------------------------------- 427 428 CALL obs_coo_grd( s ladata%nsurf, sladata%mi, sladata%mj, &429 & s ladata%nqc, igrdobs )133 ! Check for surface data failing the grid search 134 ! ----------------------------------------------------------------------- 135 136 CALL obs_coo_grd( surfdata%nsurf, surfdata%mi, surfdata%mj, & 137 & surfdata%nqc, igrdobs ) 430 138 431 139 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) … … 435 143 ! ----------------------------------------------------------------------- 436 144 437 CALL obs_coo_spc_2d( s ladata%nsurf, &145 CALL obs_coo_spc_2d( surfdata%nsurf, & 438 146 & jpi, jpj, & 439 & s ladata%mi, sladata%mj, &440 & s ladata%rlam, sladata%rphi, &147 & surfdata%mi, surfdata%mj, & 148 & surfdata%rlam, surfdata%rphi, & 441 149 & glamt, gphit, & 442 & tmask(:,:,1), s ladata%nqc, &150 & tmask(:,:,1), surfdata%nqc, & 443 151 & iosdsobs, ilansobs, & 444 152 & inlasobs, ld_nea ) … … 449 157 450 158 ! ----------------------------------------------------------------------- 451 ! Copy useful data from the s ladata data structure to452 ! the s ladatqc data structure159 ! Copy useful data from the surfdata data structure to 160 ! the surfdataqc data structure 453 161 ! ----------------------------------------------------------------------- 454 162 455 163 ! Allocate the selection arrays 456 164 457 ALLOCATE( llvalid(s ladata%nsurf) )165 ALLOCATE( llvalid(surfdata%nsurf) ) 458 166 459 167 ! We want all data which has qc flags <= 10 460 168 461 llvalid(:) = ( s ladata%nqc(:) <= 10 )169 llvalid(:) = ( surfdata%nqc(:) <= 10 ) 462 170 463 171 ! The actual copying 464 172 465 CALL obs_surf_compress( s ladata, sladatqc, .TRUE., numout, &173 CALL obs_surf_compress( surfdata, surfdataqc, .TRUE., numout, & 466 174 & lvalid=llvalid ) 467 175 … … 477 185 IF(lwp) THEN 478 186 WRITE(numout,*) 479 WRITE(numout,*) 'obs_pre_s la:'187 WRITE(numout,*) 'obs_pre_surf :' 480 188 WRITE(numout,*) '~~~~~~~~~~~' 481 189 WRITE(numout,*) 482 WRITE(numout,*) ' S LA data outside time domain= ', &190 WRITE(numout,*) ' Surface data outside time domain = ', & 483 191 & iotdobsmpp 484 WRITE(numout,*) ' Remaining SLAdata that failed grid search = ', &192 WRITE(numout,*) ' Remaining surface data that failed grid search = ', & 485 193 & igrdobsmpp 486 WRITE(numout,*) ' Remaining SLAdata outside space domain = ', &194 WRITE(numout,*) ' Remaining surface data outside space domain = ', & 487 195 & iosdsobsmpp 488 WRITE(numout,*) ' Remaining SLAdata at land points = ', &196 WRITE(numout,*) ' Remaining surface data at land points = ', & 489 197 & ilansobsmpp 490 198 IF (ld_nea) THEN 491 WRITE(numout,*) ' Remaining SLAdata near land points (removed) = ', &199 WRITE(numout,*) ' Remaining surface data near land points (removed) = ', & 492 200 & inlasobsmpp 493 201 ELSE 494 WRITE(numout,*) ' Remaining SLAdata near land points (kept) = ', &202 WRITE(numout,*) ' Remaining surface data near land points (kept) = ', & 495 203 & inlasobsmpp 496 204 ENDIF 497 WRITE(numout,*) ' SLAdata accepted = ', &498 & s ladatqc%nsurfmpp205 WRITE(numout,*) ' surface data accepted = ', & 206 & surfdataqc%nsurfmpp 499 207 500 208 WRITE(numout,*) … … 505 213 ENDIF 506 214 507 DO jobs = 1, s ladatqc%nsurf508 inrc = s ladatqc%mstp(jobs) + 2 - nit000509 s ladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1215 DO jobs = 1, surfdataqc%nsurf 216 inrc = surfdataqc%mstp(jobs) + 2 - nit000 217 surfdataqc%nsstp(inrc) = surfdataqc%nsstp(inrc) + 1 510 218 END DO 511 219 512 CALL obs_mpp_sum_integers( s ladatqc%nsstp, sladatqc%nsstpmpp, &220 CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, & 513 221 & nitend - nit000 + 2 ) 514 222 … … 516 224 DO jstp = nit000 - 1, nitend 517 225 inrc = jstp - nit000 + 2 518 WRITE(numout,1999) jstp, s ladatqc%nsstpmpp(inrc)226 WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc) 519 227 END DO 520 228 ENDIF … … 524 232 1999 FORMAT(10X,I9,5X,I17) 525 233 526 END SUBROUTINE obs_pre_sla 527 528 SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea ) 529 !!---------------------------------------------------------------------- 530 !! *** ROUTINE obs_pre_sst *** 531 !! 532 !! ** Purpose : First level check and screening of SST observations 533 !! 534 !! ** Method : First level check and screening of SST observations 535 !! 536 !! ** Action : 537 !! 538 !! References : 539 !! 540 !! History : 541 !! ! 2007-03 (S. Ricci) SST data preparation 542 !!---------------------------------------------------------------------- 543 !! * Modules used 544 USE domstp ! Domain: set the time-step 545 USE par_oce ! Ocean parameters 546 USE dom_oce, ONLY : & ! Geographical information 547 & glamt, & 548 & gphit, & 549 & tmask, & 550 & nproc 551 !! * Arguments 552 TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data 553 TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening 554 LOGICAL :: ld_sst ! Switch for SST data 555 LOGICAL :: ld_nea ! Switch for rejecting observation near land 556 !! * Local declarations 557 INTEGER :: iyea0 ! Initial date 558 INTEGER :: imon0 ! - (year, month, day, hour, minute) 559 INTEGER :: iday0 560 INTEGER :: ihou0 561 INTEGER :: imin0 562 INTEGER :: icycle ! Current assimilation cycle 563 ! Counters for observations that 564 INTEGER :: iotdobs ! - outside time domain 565 INTEGER :: iosdsobs ! - outside space domain 566 INTEGER :: ilansobs ! - within a model land cell 567 INTEGER :: inlasobs ! - close to land 568 INTEGER :: igrdobs ! - fail the grid search 569 ! Global counters for observations that 570 INTEGER :: iotdobsmpp ! - outside time domain 571 INTEGER :: iosdsobsmpp ! - outside space domain 572 INTEGER :: ilansobsmpp ! - within a model land cell 573 INTEGER :: inlasobsmpp ! - close to land 574 INTEGER :: igrdobsmpp ! - fail the grid search 575 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 576 & llvalid ! SST data selection 577 INTEGER :: jobs ! Obs. loop variable 578 INTEGER :: jstp ! Time loop variable 579 INTEGER :: inrc ! Time index variable 580 581 IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...' 582 583 ! Initial date initialization (year, month, day, hour, minute) 584 iyea0 = ndate0 / 10000 585 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 586 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 587 ihou0 = 0 588 imin0 = 0 589 590 icycle = no ! Assimilation cycle 591 592 ! Diagnotics counters for various failures. 593 594 iotdobs = 0 595 igrdobs = 0 596 iosdsobs = 0 597 ilansobs = 0 598 inlasobs = 0 599 600 ! ----------------------------------------------------------------------- 601 ! Find time coordinate for SST data 602 ! ----------------------------------------------------------------------- 603 604 CALL obs_coo_tim( icycle, & 605 & iyea0, imon0, iday0, ihou0, imin0, & 606 & sstdata%nsurf, sstdata%nyea, sstdata%nmon, & 607 & sstdata%nday, sstdata%nhou, sstdata%nmin, & 608 & sstdata%nqc, sstdata%mstp, iotdobs ) 609 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 610 ! ----------------------------------------------------------------------- 611 ! Check for SST data failing the grid search 612 ! ----------------------------------------------------------------------- 613 614 CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, & 615 & sstdata%nqc, igrdobs ) 616 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 617 618 ! ----------------------------------------------------------------------- 619 ! Check for land points. 620 ! ----------------------------------------------------------------------- 621 622 CALL obs_coo_spc_2d( sstdata%nsurf, & 623 & jpi, jpj, & 624 & sstdata%mi, sstdata%mj, & 625 & sstdata%rlam, sstdata%rphi, & 626 & glamt, gphit, & 627 & tmask(:,:,1), sstdata%nqc, & 628 & iosdsobs, ilansobs, & 629 & inlasobs, ld_nea ) 630 631 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 632 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 633 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 634 635 ! ----------------------------------------------------------------------- 636 ! Copy useful data from the sstdata data structure to 637 ! the sstdatqc data structure 638 ! ----------------------------------------------------------------------- 639 640 ! Allocate the selection arrays 641 642 ALLOCATE( llvalid(sstdata%nsurf) ) 643 644 ! We want all data which has qc flags <= 0 645 646 llvalid(:) = ( sstdata%nqc(:) <= 10 ) 647 648 ! The actual copying 649 650 CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, & 651 & lvalid=llvalid ) 652 653 ! Dellocate the selection arrays 654 DEALLOCATE( llvalid ) 655 656 ! ----------------------------------------------------------------------- 657 ! Print information about what observations are left after qc 658 ! ----------------------------------------------------------------------- 659 660 ! Update the total observation counter array 661 662 IF(lwp) THEN 663 WRITE(numout,*) 664 WRITE(numout,*) 'obs_pre_sst :' 665 WRITE(numout,*) '~~~~~~~~~~~' 666 WRITE(numout,*) 667 WRITE(numout,*) ' SST data outside time domain = ', & 668 & iotdobsmpp 669 WRITE(numout,*) ' Remaining SST data that failed grid search = ', & 670 & igrdobsmpp 671 WRITE(numout,*) ' Remaining SST data outside space domain = ', & 672 & iosdsobsmpp 673 WRITE(numout,*) ' Remaining SST data at land points = ', & 674 & ilansobsmpp 675 IF (ld_nea) THEN 676 WRITE(numout,*) ' Remaining SST data near land points (removed) = ', & 677 & inlasobsmpp 678 ELSE 679 WRITE(numout,*) ' Remaining SST data near land points (kept) = ', & 680 & inlasobsmpp 681 ENDIF 682 WRITE(numout,*) ' SST data accepted = ', & 683 & sstdatqc%nsurfmpp 684 685 WRITE(numout,*) 686 WRITE(numout,*) ' Number of observations per time step :' 687 WRITE(numout,*) 688 WRITE(numout,1997) 689 WRITE(numout,1998) 690 ENDIF 691 692 DO jobs = 1, sstdatqc%nsurf 693 inrc = sstdatqc%mstp(jobs) + 2 - nit000 694 sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1 695 END DO 696 697 CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, & 698 & nitend - nit000 + 2 ) 699 700 IF ( lwp ) THEN 701 DO jstp = nit000 - 1, nitend 702 inrc = jstp - nit000 + 2 703 WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc) 704 END DO 705 ENDIF 706 707 1997 FORMAT(10X,'Time step',5X,'Sea surface temperature') 708 1998 FORMAT(10X,'---------',5X,'-----------------') 709 1999 FORMAT(10X,I9,5X,I17) 710 711 END SUBROUTINE obs_pre_sst 712 713 SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea ) 714 !!---------------------------------------------------------------------- 715 !! *** ROUTINE obs_pre_seaice *** 716 !! 717 !! ** Purpose : First level check and screening of Sea Ice observations 718 !! 719 !! ** Method : First level check and screening of Sea Ice observations 720 !! 721 !! ** Action : 722 !! 723 !! References : 724 !! 725 !! History : 726 !! ! 2007-11 (D. Lea) based on obs_pre_sst 727 !!---------------------------------------------------------------------- 728 !! * Modules used 729 USE domstp ! Domain: set the time-step 730 USE par_oce ! Ocean parameters 731 USE dom_oce, ONLY : & ! Geographical information 732 & glamt, & 733 & gphit, & 734 & tmask, & 735 & nproc 736 !! * Arguments 737 TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data 738 TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening 739 LOGICAL :: ld_seaice ! Switch for sea ice data 740 LOGICAL :: ld_nea ! Switch for rejecting observation near land 741 !! * Local declarations 742 INTEGER :: iyea0 ! Initial date 743 INTEGER :: imon0 ! - (year, month, day, hour, minute) 744 INTEGER :: iday0 745 INTEGER :: ihou0 746 INTEGER :: imin0 747 INTEGER :: icycle ! Current assimilation cycle 748 ! Counters for observations that 749 INTEGER :: iotdobs ! - outside time domain 750 INTEGER :: iosdsobs ! - outside space domain 751 INTEGER :: ilansobs ! - within a model land cell 752 INTEGER :: inlasobs ! - close to land 753 INTEGER :: igrdobs ! - fail the grid search 754 ! Global counters for observations that 755 INTEGER :: iotdobsmpp ! - outside time domain 756 INTEGER :: iosdsobsmpp ! - outside space domain 757 INTEGER :: ilansobsmpp ! - within a model land cell 758 INTEGER :: inlasobsmpp ! - close to land 759 INTEGER :: igrdobsmpp ! - fail the grid search 760 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 761 & llvalid ! data selection 762 INTEGER :: jobs ! Obs. loop variable 763 INTEGER :: jstp ! Time loop variable 764 INTEGER :: inrc ! Time index variable 765 766 IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...' 767 768 ! Initial date initialization (year, month, day, hour, minute) 769 iyea0 = ndate0 / 10000 770 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 771 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 772 ihou0 = 0 773 imin0 = 0 774 775 icycle = no ! Assimilation cycle 776 777 ! Diagnotics counters for various failures. 778 779 iotdobs = 0 780 igrdobs = 0 781 iosdsobs = 0 782 ilansobs = 0 783 inlasobs = 0 784 785 ! ----------------------------------------------------------------------- 786 ! Find time coordinate for sea ice data 787 ! ----------------------------------------------------------------------- 788 789 CALL obs_coo_tim( icycle, & 790 & iyea0, imon0, iday0, ihou0, imin0, & 791 & seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, & 792 & seaicedata%nday, seaicedata%nhou, seaicedata%nmin, & 793 & seaicedata%nqc, seaicedata%mstp, iotdobs ) 794 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 795 ! ----------------------------------------------------------------------- 796 ! Check for sea ice data failing the grid search 797 ! ----------------------------------------------------------------------- 798 799 CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, & 800 & seaicedata%nqc, igrdobs ) 801 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 802 803 ! ----------------------------------------------------------------------- 804 ! Check for land points. 805 ! ----------------------------------------------------------------------- 806 807 CALL obs_coo_spc_2d( seaicedata%nsurf, & 808 & jpi, jpj, & 809 & seaicedata%mi, seaicedata%mj, & 810 & seaicedata%rlam, seaicedata%rphi, & 811 & glamt, gphit, & 812 & tmask(:,:,1), seaicedata%nqc, & 813 & iosdsobs, ilansobs, & 814 & inlasobs, ld_nea ) 815 816 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 817 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 818 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 819 820 ! ----------------------------------------------------------------------- 821 ! Copy useful data from the seaicedata data structure to 822 ! the seaicedatqc data structure 823 ! ----------------------------------------------------------------------- 824 825 ! Allocate the selection arrays 826 827 ALLOCATE( llvalid(seaicedata%nsurf) ) 828 829 ! We want all data which has qc flags <= 0 830 831 llvalid(:) = ( seaicedata%nqc(:) <= 10 ) 832 833 ! The actual copying 834 835 CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, & 836 & lvalid=llvalid ) 837 838 ! Dellocate the selection arrays 839 DEALLOCATE( llvalid ) 840 841 ! ----------------------------------------------------------------------- 842 ! Print information about what observations are left after qc 843 ! ----------------------------------------------------------------------- 844 845 ! Update the total observation counter array 846 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'obs_pre_seaice :' 850 WRITE(numout,*) '~~~~~~~~~~~' 851 WRITE(numout,*) 852 WRITE(numout,*) ' Sea ice data outside time domain = ', & 853 & iotdobsmpp 854 WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', & 855 & igrdobsmpp 856 WRITE(numout,*) ' Remaining sea ice data outside space domain = ', & 857 & iosdsobsmpp 858 WRITE(numout,*) ' Remaining sea ice data at land points = ', & 859 & ilansobsmpp 860 IF (ld_nea) THEN 861 WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', & 862 & inlasobsmpp 863 ELSE 864 WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', & 865 & inlasobsmpp 866 ENDIF 867 WRITE(numout,*) ' Sea ice data accepted = ', & 868 & seaicedatqc%nsurfmpp 869 870 WRITE(numout,*) 871 WRITE(numout,*) ' Number of observations per time step :' 872 WRITE(numout,*) 873 WRITE(numout,1997) 874 WRITE(numout,1998) 875 ENDIF 876 877 DO jobs = 1, seaicedatqc%nsurf 878 inrc = seaicedatqc%mstp(jobs) + 2 - nit000 879 seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1 880 END DO 881 882 CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, & 883 & nitend - nit000 + 2 ) 884 885 IF ( lwp ) THEN 886 DO jstp = nit000 - 1, nitend 887 inrc = jstp - nit000 + 2 888 WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc) 889 END DO 890 ENDIF 891 892 1997 FORMAT(10X,'Time step',5X,'Sea ice data ') 893 1998 FORMAT(10X,'---------',5X,'-----------------') 894 1999 FORMAT(10X,I9,5X,I17) 895 896 END SUBROUTINE obs_pre_seaice 897 898 SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav ) 899 !!---------------------------------------------------------------------- 900 !! *** ROUTINE obs_pre_taovel *** 901 !! 902 !! ** Purpose : First level check and screening of U and V profiles 903 !! 904 !! ** Method : First level check and screening of U and V profiles 234 END SUBROUTINE obs_pre_surf 235 236 237 SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav ) 238 !!---------------------------------------------------------------------- 239 !! *** ROUTINE obs_pre_prof *** 240 !! 241 !! ** Purpose : First level check and screening of profiles 242 !! 243 !! ** Method : First level check and screening of profiles 905 244 !! 906 245 !! History : … … 908 247 !! ! 2008-09 (M. Valdivieso) : TAO velocity data 909 248 !! ! 2009-01 (K. Mogensen) : New feedback strictuer 249 !! ! 2015-02 (M. Martin) : Combined profile routine. 910 250 !! 911 251 !!---------------------------------------------------------------------- … … 962 302 INTEGER :: inrc ! Time index variable 963 303 964 IF(lwp) WRITE(numout,*)'obs_pre_ vel: Preparing the velocityprofile data'304 IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data' 965 305 966 306 ! Initial date initialization (year, month, day, hour, minute) … … 1186 526 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 1187 527 1188 END SUBROUTINE obs_pre_ vel528 END SUBROUTINE obs_pre_prof 1189 529 1190 530 SUBROUTINE obs_coo_tim( kcycle, & -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r4990 r5659 33 33 PRIVATE 34 34 35 PUBLIC obs_rea_pro _dri! Read the profile observations35 PUBLIC obs_rea_prof ! Read the profile observations 36 36 37 37 !!---------------------------------------------------------------------- … … 43 43 CONTAINS 44 44 45 SUBROUTINE obs_rea_pro_dri( kformat, & 46 & profdata, knumfiles, cfilenames, & 47 & kvars, kextr, kstp, ddobsini, ddobsend, & 48 & ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, & 49 & ldmod, kdailyavtypes ) 45 SUBROUTINE obs_rea_prof( profdata, knumfiles, cfilenames, & 46 & kvars, kextr, kstp, ddobsini, ddobsend, & 47 & ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, & 48 & ldmod, kdailyavtypes ) 50 49 !!--------------------------------------------------------------------- 51 50 !! 52 !! *** ROUTINE obs_rea_pro _dri***51 !! *** ROUTINE obs_rea_prof *** 53 52 !! 54 53 !! ** Purpose : Read from file the profile observations 55 54 !! 56 !! ** Method : Depending on kformat either ENACT, CORIOLIS or57 !! feedback data files are read55 !! ** Method : Read feedback data in and transform to NEMO internal 56 !! profile data structure 58 57 !! 59 58 !! ** Action : … … 67 66 68 67 !! * Arguments 69 INTEGER :: kformat ! Format of input data70 ! ! 1: ENACT71 ! ! 2: Coriolis72 68 TYPE(obs_prof), INTENT(OUT) :: profdata ! Profile data to be read 73 69 INTEGER, INTENT(IN) :: knumfiles ! Number of files to read in … … 89 85 90 86 !! * Local declarations 91 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_pro _dri'87 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 92 88 INTEGER :: jvar 93 89 INTEGER :: ji … … 195 191 !--------------------------------------------------------------------- 196 192 197 iflag = nf90_open( TRIM( TRIM( cfilenames(jj)) ), nf90_nowrite, &193 iflag = nf90_open( TRIM( cfilenames(jj) ), nf90_nowrite, & 198 194 & i_file_id ) 199 195 … … 202 198 IF ( ldignmis ) THEN 203 199 inpfiles(jj)%nobs = 0 204 CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj)) ) // &200 CALL ctl_warn( 'File ' // TRIM( cfilenames(jj) ) // & 205 201 & ' not found' ) 206 202 ELSE 207 CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj)) ) // &203 CALL ctl_stop( 'File ' // TRIM( cfilenames(jj) ) // & 208 204 & ' not found' ) 209 205 ENDIF … … 220 216 ! Read the profile file into inpfiles 221 217 !------------------------------------------------------------------ 222 IF ( kformat == 0 ) THEN 223 CALL init_obfbdata( inpfiles(jj) ) 224 IF(lwp) THEN 225 WRITE(numout,*) 226 WRITE(numout,*)'Reading from feedback file :', & 227 & TRIM( cfilenames(jj) ) 228 ENDIF 229 CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), & 230 & ldgrid = .TRUE. ) 231 IF ( inpfiles(jj)%nvar < 2 ) THEN 232 CALL ctl_stop( 'Feedback format error' ) 233 RETURN 234 ENDIF 235 IF ( TRIM(inpfiles(jj)%cname(1)) /= 'POTM' ) THEN 236 CALL ctl_stop( 'Feedback format error' ) 237 RETURN 238 ENDIF 239 IF ( TRIM(inpfiles(jj)%cname(2)) /= 'PSAL' ) THEN 240 CALL ctl_stop( 'Feedback format error' ) 241 RETURN 242 ENDIF 243 IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 244 CALL ctl_stop( 'Model not in input data' ) 245 RETURN 246 ENDIF 247 ELSEIF ( kformat == 1 ) THEN 248 CALL read_enactfile( TRIM( cfilenames(jj) ), inpfiles(jj), & 249 & numout, lwp, .TRUE. ) 250 ELSEIF ( kformat == 2 ) THEN 251 CALL read_coriofile( TRIM( cfilenames(jj) ), inpfiles(jj), & 252 & numout, lwp, .TRUE. ) 253 ELSE 254 CALL ctl_stop( 'File format unknown' ) 218 CALL init_obfbdata( inpfiles(jj) ) 219 IF(lwp) THEN 220 WRITE(numout,*) 221 WRITE(numout,*)'Reading from feedback file :', & 222 & TRIM( cfilenames(jj) ) 223 ENDIF 224 CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), & 225 & ldgrid = .TRUE. ) 226 227 IF ( inpfiles(jj)%nvar < 2 ) THEN 228 CALL ctl_stop( 'Feedback format error' ) 229 RETURN 230 ENDIF 231 IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 232 CALL ctl_stop( 'Model not in input data' ) 233 RETURN 255 234 ENDIF 256 235 … … 836 815 DEALLOCATE( inpfiles ) 837 816 838 END SUBROUTINE obs_rea_pro _dri817 END SUBROUTINE obs_rea_prof 839 818 840 819 END MODULE obs_read_prof -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
r4990 r5659 6 6 7 7 !!---------------------------------------------------------------------- 8 !! obs_wri_p3d : Write profile observation diagnostics in NetCDF format 9 !! obs_wri_sla : Write SLA observation related diagnostics 10 !! obs_wri_sst : Write SST observation related diagnostics 11 !! obs_wri_seaice: Write seaice observation related diagnostics 12 !! obs_wri_vel : Write velocity observation diagnostics in NetCDF format 8 !! obs_wri_prof : Write profile observations in feedback format 9 !! obs_wri_surf : Write surface observations in feedback format 13 10 !! obs_wri_stats : Print basic statistics on the data being written out 14 11 !!---------------------------------------------------------------------- … … 30 27 USE obs_conv ! Conversion between units 31 28 USE obs_const 32 USE obs_sla_types33 29 USE obs_rot_vel ! Rotation of velocities 34 30 USE obs_mpp ! MPP support routines for observation diagnostics … … 39 35 !! * Routine accessibility 40 36 PRIVATE 41 PUBLIC obs_wri_p3d, & ! Write profile observation related diagnostics 42 & obs_wri_sla, & ! Write SLA observation related diagnostics 43 & obs_wri_sst, & ! Write SST observation related diagnostics 44 & obs_wri_sss, & ! Write SSS observation related diagnostics 45 & obs_wri_seaice, & ! Write seaice observation related diagnostics 46 & obs_wri_vel, & ! Write velocity observation related diagnostics 37 PUBLIC obs_wri_prof, & ! Write profile observation files 38 & obs_wri_surf, & ! Write surface observation files 47 39 & obswriinfo 48 40 … … 63 55 CONTAINS 64 56 65 SUBROUTINE obs_wri_p 3d( cprefix, profdata, padd, pext )57 SUBROUTINE obs_wri_prof( cobstype, profdata, padd, pext ) 66 58 !!----------------------------------------------------------------------- 67 59 !! 68 !! *** ROUTINE obs_wri_p3d *** 69 !! 70 !! ** Purpose : Write temperature and salinity (profile) observation 71 !! related diagnostics 60 !! *** ROUTINE obs_wri_prof *** 61 !! 62 !! ** Purpose : Write profile feedback files 72 63 !! 73 64 !! ** Method : NetCDF … … 82 73 !! ! 07-03 (K. Mogensen) General handling of profiles 83 74 !! ! 09-01 (K. Mogensen) New feedback format 75 !! ! 15-02 (M. Martin) Combined routine for writing profiles 84 76 !!----------------------------------------------------------------------- 85 77 … … 87 79 88 80 !! * Arguments 89 CHARACTER(LEN=*), INTENT(IN) :: c prefix! Prefix for output files81 CHARACTER(LEN=*), INTENT(IN) :: cobstype ! Prefix for output files 90 82 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 91 83 TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable … … 125 117 ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 126 118 END DO 127 CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 128 & 1 + nadd, 1 + next, .TRUE. ) 129 130 fbdata%cname(1) = 'POTM' 131 fbdata%cname(2) = 'PSAL' 132 fbdata%coblong(1) = 'Potential temperature' 133 fbdata%coblong(2) = 'Practical salinity' 134 fbdata%cobunit(1) = 'Degrees centigrade' 135 fbdata%cobunit(2) = 'PSU' 136 fbdata%cextname(1) = 'TEMP' 137 fbdata%cextlong(1) = 'Insitu temperature' 138 fbdata%cextunit(1) = 'Degrees centigrade' 139 DO je = 1, next 140 fbdata%cextname(1+je) = pext%cdname(je) 141 fbdata%cextlong(1+je) = pext%cdlong(je,1) 142 fbdata%cextunit(1+je) = pext%cdunit(je,1) 143 END DO 119 120 SELECT CASE ( TRIM(cobstype) ) 121 CASE('prof') 122 123 CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 124 & 1 + nadd, 1 + next, .TRUE. ) 125 fbdata%cname(1) = 'POTM' 126 fbdata%cname(2) = 'PSAL' 127 fbdata%coblong(1) = 'Potential temperature' 128 fbdata%coblong(2) = 'Practical salinity' 129 fbdata%cobunit(1) = 'Degrees centigrade' 130 fbdata%cobunit(2) = 'PSU' 131 fbdata%cextname(1) = 'TEMP' 132 fbdata%cextlong(1) = 'Insitu temperature' 133 fbdata%cextunit(1) = 'Degrees centigrade' 134 fbdata%caddlong(1,1) = 'Model interpolated potential temperature' 135 fbdata%caddlong(1,2) = 'Model interpolated practical salinity' 136 fbdata%caddunit(1,1) = 'Degrees centigrade' 137 fbdata%caddunit(1,2) = 'PSU' 138 fbdata%cgrid(:) = 'T' 139 DO je = 1, next 140 fbdata%cextname(1+je) = pext%cdname(je) 141 fbdata%cextlong(1+je) = pext%cdlong(je,1) 142 fbdata%cextunit(1+je) = pext%cdunit(je,1) 143 END DO 144 DO ja = 1, nadd 145 fbdata%caddname(1+ja) = padd%cdname(ja) 146 DO jvar = 1, 2 147 fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar) 148 fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar) 149 END DO 150 END DO 151 152 CASE('vel') 153 154 CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. ) 155 fbdata%cname(1) = 'UVEL' 156 fbdata%cname(2) = 'VVEL' 157 fbdata%coblong(1) = 'Zonal velocity' 158 fbdata%coblong(2) = 'Meridional velocity' 159 fbdata%cobunit(1) = 'm/s' 160 fbdata%cobunit(2) = 'm/s' 161 DO je = 1, next 162 fbdata%cextname(je) = pext%cdname(je) 163 fbdata%cextlong(je) = pext%cdlong(je,1) 164 fbdata%cextunit(je) = pext%cdunit(je,1) 165 END DO 166 fbdata%caddlong(1,1) = 'Model interpolated zonal velocity' 167 fbdata%caddlong(1,2) = 'Model interpolated meridional velocity' 168 fbdata%caddunit(1,1) = 'm/s' 169 fbdata%caddunit(1,2) = 'm/s' 170 fbdata%caddname(2) = 'HxG' 171 fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)' 172 fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)' 173 fbdata%caddunit(2,1) = 'm/s' 174 fbdata%caddunit(2,2) = 'm/s' 175 fbdata%cgrid(1) = 'U' 176 fbdata%cgrid(2) = 'V' 177 DO ja = 1, nadd 178 fbdata%caddname(2+ja) = padd%cdname(ja) 179 fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 180 fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 181 END DO 182 ALLOCATE( & 183 & zu(profdata%nvprot(1)), & 184 & zv(profdata%nvprot(2)) & 185 & ) 186 CALL obs_rotvel( profdata, k2dint, zu, zv ) 187 188 END SELECT 189 144 190 fbdata%caddname(1) = 'Hx' 145 fbdata%caddlong(1,1) = 'Model interpolated potential temperature'146 fbdata%caddlong(1,2) = 'Model interpolated practical salinity'147 fbdata%caddunit(1,1) = 'Degrees centigrade'148 fbdata%caddunit(1,2) = 'PSU'149 fbdata%cgrid(:) = 'T'150 DO ja = 1, nadd151 fbdata%caddname(1+ja) = padd%cdname(ja)152 DO jvar = 1, 2153 fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)154 fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)155 END DO156 END DO157 191 158 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(c prefix), nproc192 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc 159 193 160 194 IF(lwp) THEN 161 195 WRITE(numout,*) 162 WRITE(numout,*)'obs_wri_p 3d:'196 WRITE(numout,*)'obs_wri_prof :' 163 197 WRITE(numout,*)'~~~~~~~~~~~~~' 164 198 WRITE(numout,*)'Writing profile feedback file : ',TRIM(cfname) … … 208 242 DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 209 243 ik = profdata%var(jvar)%nvlidx(jk) 210 fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 244 IF ( TRIM(cobstype) == 'prof' ) THEN 245 fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 246 ELSE IF ( TRIM(cobstype) == 'vel' ) THEN 247 IF ( jvar == 1 ) THEN 248 fbdata%padd(ik,jo,1,jvar) = zu(jk) 249 ELSE 250 fbdata%padd(ik,jo,1,jvar) = zv(jk) 251 ENDIF 252 fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk) 253 ENDIF 211 254 fbdata%pob(ik,jo,jvar) = profdata%var(jvar)%vobs(jk) 212 255 fbdata%pdep(ik,jo) = profdata%var(jvar)%vdep(jk) … … 237 280 END DO 238 281 239 ! Convert insitu temperature to potential temperature using the model 240 ! salinity if no potential temperature 241 DO jo = 1, fbdata%nobs 242 IF ( fbdata%pphi(jo) < 9999.0 ) THEN 243 DO jk = 1, fbdata%nlev 244 IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. & 245 & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. & 246 & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. & 247 & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN 248 zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), & 249 & REAL(fbdata%pphi(jo),wp) ) 250 fbdata%pob(jk,jo,1) = potemp( & 251 & REAL(fbdata%padd(jk,jo,1,2), wp), & 252 & REAL(fbdata%pext(jk,jo,1), wp), & 253 & zpres, 0.0_wp ) 254 ENDIF 255 END DO 256 ENDIF 257 END DO 282 IF ( TRIM(cobstype) == 'prof' ) THEN 283 ! Convert insitu temperature to potential temperature using the model 284 ! salinity if no potential temperature 285 DO jo = 1, fbdata%nobs 286 IF ( fbdata%pphi(jo) < 9999.0 ) THEN 287 DO jk = 1, fbdata%nlev 288 IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. & 289 & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. & 290 & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. & 291 & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN 292 zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), & 293 & REAL(fbdata%pphi(jo),wp) ) 294 fbdata%pob(jk,jo,1) = potemp( & 295 & REAL(fbdata%padd(jk,jo,1,2), wp), & 296 & REAL(fbdata%pext(jk,jo,1), wp), & 297 & zpres, 0.0_wp ) 298 ENDIF 299 END DO 300 ENDIF 301 END DO 302 ENDIF 258 303 259 304 ! Write the obfbdata structure … … 265 310 CALL dealloc_obfbdata( fbdata ) 266 311 267 END SUBROUTINE obs_wri_p 3d268 269 SUBROUTINE obs_wri_s la( cprefix, sladata, padd, pext )312 END SUBROUTINE obs_wri_prof 313 314 SUBROUTINE obs_wri_surf( cobstype, surfdata, padd, pext ) 270 315 !!----------------------------------------------------------------------- 271 316 !! 272 !! *** ROUTINE obs_wri_sla *** 273 !! 274 !! ** Purpose : Write SLA observation diagnostics 275 !! related 317 !! *** ROUTINE obs_wri_surf *** 318 !! 319 !! ** Purpose : Write surface observation files 276 320 !! 277 321 !! ** Method : NetCDF … … 281 325 !! ! 07-03 (K. Mogensen) Original 282 326 !! ! 09-01 (K. Mogensen) New feedback format. 327 !! ! 15-02 (M. Martin) Combined surface writing routine. 283 328 !!----------------------------------------------------------------------- 284 329 … … 287 332 288 333 !! * Arguments 289 CHARACTER(LEN=*), INTENT(IN) :: c prefix! Prefix for output files290 TYPE(obs_surf), INTENT(INOUT) :: s ladata ! Full set of SLAa334 CHARACTER(LEN=*), INTENT(IN) :: cobstype ! Prefix for output files 335 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data 291 336 TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable 292 337 TYPE(obswriinfo), OPTIONAL :: pext ! Extra info … … 295 340 TYPE(obfbdata) :: fbdata 296 341 CHARACTER(LEN=40) :: cfname ! netCDF filename 297 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_s la'342 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 298 343 INTEGER :: jo 299 344 INTEGER :: ja … … 316 361 CALL init_obfbdata( fbdata ) 317 362 318 CALL alloc_obfbdata( fbdata, 1, s ladata%nsurf, 1, &363 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 319 364 & 2 + nadd, 1 + next, .TRUE. ) 320 365 321 fbdata%cname(1) = 'SLA' 322 fbdata%coblong(1) = 'Sea level anomaly' 323 fbdata%cobunit(1) = 'Metres' 324 fbdata%cextname(1) = 'MDT' 325 fbdata%cextlong(1) = 'Mean dynamic topography' 326 fbdata%cextunit(1) = 'Metres' 327 DO je = 1, next 328 fbdata%cextname(1+je) = pext%cdname(je) 329 fbdata%cextlong(1+je) = pext%cdlong(je,1) 330 fbdata%cextunit(1+je) = pext%cdunit(je,1) 331 END DO 366 SELECT CASE ( TRIM(cobstype) ) 367 CASE('sla') 368 369 fbdata%cname(1) = 'SLA' 370 fbdata%coblong(1) = 'Sea level anomaly' 371 fbdata%cobunit(1) = 'Metres' 372 fbdata%cextname(1) = 'MDT' 373 fbdata%cextlong(1) = 'Mean dynamic topography' 374 fbdata%cextunit(1) = 'Metres' 375 DO je = 1, next 376 fbdata%cextname(je) = pext%cdname(je) 377 fbdata%cextlong(je) = pext%cdlong(je,1) 378 fbdata%cextunit(je) = pext%cdunit(je,1) 379 END DO 380 fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT' 381 fbdata%caddunit(1,1) = 'Metres' 382 fbdata%caddname(2) = 'SSH' 383 fbdata%caddlong(2,1) = 'Model Sea surface height' 384 fbdata%caddunit(2,1) = 'Metres' 385 fbdata%cgrid(1) = 'T' 386 DO ja = 1, nadd 387 fbdata%caddname(2+ja) = padd%cdname(ja) 388 fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 389 fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 390 END DO 391 392 CASE('sst') 393 394 fbdata%cname(1) = 'SST' 395 fbdata%coblong(1) = 'Sea surface temperature' 396 fbdata%cobunit(1) = 'Degree centigrade' 397 DO je = 1, next 398 fbdata%cextname(je) = pext%cdname(je) 399 fbdata%cextlong(je) = pext%cdlong(je,1) 400 fbdata%cextunit(je) = pext%cdunit(je,1) 401 END DO 402 fbdata%caddlong(1,1) = 'Model interpolated SST' 403 fbdata%caddunit(1,1) = 'Degree centigrade' 404 fbdata%cgrid(1) = 'T' 405 DO ja = 1, nadd 406 fbdata%caddname(1+ja) = padd%cdname(ja) 407 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 408 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 409 END DO 410 411 CASE('sss') 412 413 fbdata%cname(1) = 'SSS' 414 fbdata%coblong(1) = 'Sea surface salinity' 415 fbdata%cobunit(1) = 'PSU' 416 DO je = 1, next 417 fbdata%cextname(je) = pext%cdname(je) 418 fbdata%cextlong(je) = pext%cdlong(je,1) 419 fbdata%cextunit(je) = pext%cdunit(je,1) 420 END DO 421 fbdata%caddlong(1,1) = 'Model interpolated SSS' 422 fbdata%caddunit(1,1) = 'PSU' 423 fbdata%cgrid(1) = 'T' 424 DO ja = 1, nadd 425 fbdata%caddname(1+ja) = padd%cdname(ja) 426 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 427 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 428 END DO 429 430 CASE('seaice') 431 432 fbdata%cname(1) = 'SEAICE' 433 fbdata%coblong(1) = 'Sea ice' 434 fbdata%cobunit(1) = 'Fraction' 435 DO je = 1, next 436 fbdata%cextname(je) = pext%cdname(je) 437 fbdata%cextlong(je) = pext%cdlong(je,1) 438 fbdata%cextunit(je) = pext%cdunit(je,1) 439 END DO 440 fbdata%caddlong(1,1) = 'Model interpolated ICE' 441 fbdata%caddunit(1,1) = 'Fraction' 442 fbdata%cgrid(1) = 'T' 443 DO ja = 1, nadd 444 fbdata%caddname(1+ja) = padd%cdname(ja) 445 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 446 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 447 END DO 448 449 END SELECT 450 332 451 fbdata%caddname(1) = 'Hx' 333 fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT' 334 fbdata%caddunit(1,1) = 'Metres' 335 fbdata%caddname(2) = 'SSH' 336 fbdata%caddlong(2,1) = 'Model Sea surface height' 337 fbdata%caddunit(2,1) = 'Metres' 338 fbdata%cgrid(1) = 'T' 339 DO ja = 1, nadd 340 fbdata%caddname(2+ja) = padd%cdname(ja) 341 fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 342 fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 343 END DO 344 345 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 452 453 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc 346 454 347 455 IF(lwp) THEN 348 456 WRITE(numout,*) 349 WRITE(numout,*)'obs_wri_s la:'457 WRITE(numout,*)'obs_wri_surf :' 350 458 WRITE(numout,*)'~~~~~~~~~~~~~' 351 WRITE(numout,*)'Writing SLAfeedback file : ',TRIM(cfname)459 WRITE(numout,*)'Writing surface feedback file : ',TRIM(cfname) 352 460 ENDIF 353 461 354 462 ! Transform obs_prof data structure into obfbdata structure 355 463 fbdata%cdjuldref = '19500101000000' 356 DO jo = 1, s ladata%nsurf357 fbdata%plam(jo) = s ladata%rlam(jo)358 fbdata%pphi(jo) = s ladata%rphi(jo)359 WRITE(fbdata%cdtyp(jo),'(I4)') s ladata%ntyp(jo)464 DO jo = 1, surfdata%nsurf 465 fbdata%plam(jo) = surfdata%rlam(jo) 466 fbdata%pphi(jo) = surfdata%rphi(jo) 467 WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo) 360 468 fbdata%ivqc(jo,:) = 0 361 469 fbdata%ivqcf(:,jo,:) = 0 362 IF ( s ladata%nqc(jo) > 10 ) THEN470 IF ( surfdata%nqc(jo) > 10 ) THEN 363 471 fbdata%ioqc(jo) = 4 364 472 fbdata%ioqcf(1,jo) = 0 365 fbdata%ioqcf(2,jo) = s ladata%nqc(jo) - 10473 fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 366 474 ELSE 367 fbdata%ioqc(jo) = s ladata%nqc(jo)475 fbdata%ioqc(jo) = surfdata%nqc(jo) 368 476 fbdata%ioqcf(:,jo) = 0 369 477 ENDIF … … 372 480 fbdata%itqc(jo) = 0 373 481 fbdata%itqcf(:,jo) = 0 374 fbdata%cdwmo(jo) = s ladata%cwmo(jo)375 fbdata%kindex(jo) = s ladata%nsfil(jo)482 fbdata%cdwmo(jo) = surfdata%cwmo(jo) 483 fbdata%kindex(jo) = surfdata%nsfil(jo) 376 484 IF (ln_grid_global) THEN 377 fbdata%iobsi(jo,1) = s ladata%mi(jo)378 fbdata%iobsj(jo,1) = s ladata%mj(jo)485 fbdata%iobsi(jo,1) = surfdata%mi(jo) 486 fbdata%iobsj(jo,1) = surfdata%mj(jo) 379 487 ELSE 380 fbdata%iobsi(jo,1) = mig(s ladata%mi(jo))381 fbdata%iobsj(jo,1) = mjg(s ladata%mj(jo))488 fbdata%iobsi(jo,1) = mig(surfdata%mi(jo)) 489 fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo)) 382 490 ENDIF 383 491 CALL greg2jul( 0, & 384 & s ladata%nmin(jo), &385 & s ladata%nhou(jo), &386 & s ladata%nday(jo), &387 & s ladata%nmon(jo), &388 & s ladata%nyea(jo), &492 & surfdata%nmin(jo), & 493 & surfdata%nhou(jo), & 494 & surfdata%nday(jo), & 495 & surfdata%nmon(jo), & 496 & surfdata%nyea(jo), & 389 497 & fbdata%ptim(jo), & 390 498 & krefdate = 19500101 ) 391 fbdata%padd(1,jo,1,1) = s ladata%rmod(jo,1)392 fbdata%padd(1,jo,2,1) = sladata%rext(jo,1)393 fbdata%pob(1,jo,1) = s ladata%robs(jo,1)499 fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1) 500 IF ( TRIM(cobstype) == 'sla' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1) 501 fbdata%pob(1,jo,1) = surfdata%robs(jo,1) 394 502 fbdata%pdep(1,jo) = 0.0 395 503 fbdata%idqc(1,jo) = 0 396 504 fbdata%idqcf(:,1,jo) = 0 397 IF ( s ladata%nqc(jo) > 10 ) THEN505 IF ( surfdata%nqc(jo) > 10 ) THEN 398 506 fbdata%ivqc(jo,1) = 4 399 507 fbdata%ivlqc(1,jo,1) = 4 400 508 fbdata%ivlqcf(1,1,jo,1) = 0 401 fbdata%ivlqcf(2,1,jo,1) = s ladata%nqc(jo) - 10509 fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 402 510 ELSE 403 fbdata%ivqc(jo,1) = s ladata%nqc(jo)404 fbdata%ivlqc(1,jo,1) = s ladata%nqc(jo)511 fbdata%ivqc(jo,1) = surfdata%nqc(jo) 512 fbdata%ivlqc(1,jo,1) = surfdata%nqc(jo) 405 513 fbdata%ivlqcf(:,1,jo,1) = 0 406 514 ENDIF 407 515 fbdata%iobsk(1,jo,1) = 0 408 fbdata%pext(1,jo,1) = sladata%rext(jo,2)516 IF ( TRIM(cobstype) == 'sla' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 409 517 DO ja = 1, nadd 410 518 fbdata%padd(1,jo,2+ja,1) = & 411 & s ladata%rext(jo,padd%ipoint(ja))519 & surfdata%rext(jo,padd%ipoint(ja)) 412 520 END DO 413 521 DO je = 1, next 414 522 fbdata%pext(1,jo,1+je) = & 415 & s ladata%rext(jo,pext%ipoint(je))523 & surfdata%rext(jo,pext%ipoint(je)) 416 524 END DO 417 525 END DO … … 425 533 CALL dealloc_obfbdata( fbdata ) 426 534 427 END SUBROUTINE obs_wri_sla 428 429 SUBROUTINE obs_wri_sst( cprefix, sstdata, padd, pext ) 430 !!----------------------------------------------------------------------- 431 !! 432 !! *** ROUTINE obs_wri_sst *** 433 !! 434 !! ** Purpose : Write SST observation diagnostics 435 !! related 436 !! 437 !! ** Method : NetCDF 438 !! 439 !! ** Action : 440 !! 441 !! ! 07-07 (S. Ricci) Original 442 !! ! 09-01 (K. Mogensen) New feedback format. 443 !!----------------------------------------------------------------------- 444 445 !! * Modules used 446 IMPLICIT NONE 447 448 !! * Arguments 449 CHARACTER(LEN=*), INTENT(IN) :: cprefix ! Prefix for output files 450 TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST 451 TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable 452 TYPE(obswriinfo), OPTIONAL :: pext ! Extra info 453 454 !! * Local declarations 455 TYPE(obfbdata) :: fbdata 456 CHARACTER(LEN=40) :: cfname ! netCDF filename 457 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_sst' 458 INTEGER :: jo 459 INTEGER :: ja 460 INTEGER :: je 461 INTEGER :: nadd 462 INTEGER :: next 463 464 IF ( PRESENT( padd ) ) THEN 465 nadd = padd%inum 466 ELSE 467 nadd = 0 468 ENDIF 469 470 IF ( PRESENT( pext ) ) THEN 471 next = pext%inum 472 ELSE 473 next = 0 474 ENDIF 475 476 CALL init_obfbdata( fbdata ) 477 478 CALL alloc_obfbdata( fbdata, 1, sstdata%nsurf, 1, & 479 & 1 + nadd, next, .TRUE. ) 480 481 fbdata%cname(1) = 'SST' 482 fbdata%coblong(1) = 'Sea surface temperature' 483 fbdata%cobunit(1) = 'Degree centigrade' 484 DO je = 1, next 485 fbdata%cextname(je) = pext%cdname(je) 486 fbdata%cextlong(je) = pext%cdlong(je,1) 487 fbdata%cextunit(je) = pext%cdunit(je,1) 488 END DO 489 fbdata%caddname(1) = 'Hx' 490 fbdata%caddlong(1,1) = 'Model interpolated SST' 491 fbdata%caddunit(1,1) = 'Degree centigrade' 492 fbdata%cgrid(1) = 'T' 493 DO ja = 1, nadd 494 fbdata%caddname(1+ja) = padd%cdname(ja) 495 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 496 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 497 END DO 498 499 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 500 501 IF(lwp) THEN 502 WRITE(numout,*) 503 WRITE(numout,*)'obs_wri_sst :' 504 WRITE(numout,*)'~~~~~~~~~~~~~' 505 WRITE(numout,*)'Writing SST feedback file : ',TRIM(cfname) 506 ENDIF 507 508 ! Transform obs_prof data structure into obfbdata structure 509 fbdata%cdjuldref = '19500101000000' 510 DO jo = 1, sstdata%nsurf 511 fbdata%plam(jo) = sstdata%rlam(jo) 512 fbdata%pphi(jo) = sstdata%rphi(jo) 513 WRITE(fbdata%cdtyp(jo),'(I4)') sstdata%ntyp(jo) 514 fbdata%ivqc(jo,:) = 0 515 fbdata%ivqcf(:,jo,:) = 0 516 IF ( sstdata%nqc(jo) > 10 ) THEN 517 fbdata%ioqc(jo) = 4 518 fbdata%ioqcf(1,jo) = 0 519 fbdata%ioqcf(2,jo) = sstdata%nqc(jo) - 10 520 ELSE 521 fbdata%ioqc(jo) = MAX(sstdata%nqc(jo),1) 522 fbdata%ioqcf(:,jo) = 0 523 ENDIF 524 fbdata%ipqc(jo) = 0 525 fbdata%ipqcf(:,jo) = 0 526 fbdata%itqc(jo) = 0 527 fbdata%itqcf(:,jo) = 0 528 fbdata%cdwmo(jo) = '' 529 fbdata%kindex(jo) = sstdata%nsfil(jo) 530 IF (ln_grid_global) THEN 531 fbdata%iobsi(jo,1) = sstdata%mi(jo) 532 fbdata%iobsj(jo,1) = sstdata%mj(jo) 533 ELSE 534 fbdata%iobsi(jo,1) = mig(sstdata%mi(jo)) 535 fbdata%iobsj(jo,1) = mjg(sstdata%mj(jo)) 536 ENDIF 537 CALL greg2jul( 0, & 538 & sstdata%nmin(jo), & 539 & sstdata%nhou(jo), & 540 & sstdata%nday(jo), & 541 & sstdata%nmon(jo), & 542 & sstdata%nyea(jo), & 543 & fbdata%ptim(jo), & 544 & krefdate = 19500101 ) 545 fbdata%padd(1,jo,1,1) = sstdata%rmod(jo,1) 546 fbdata%pob(1,jo,1) = sstdata%robs(jo,1) 547 fbdata%pdep(1,jo) = 0.0 548 fbdata%idqc(1,jo) = 0 549 fbdata%idqcf(:,1,jo) = 0 550 IF ( sstdata%nqc(jo) > 10 ) THEN 551 fbdata%ivqc(jo,1) = 4 552 fbdata%ivlqc(1,jo,1) = 4 553 fbdata%ivlqcf(1,1,jo,1) = 0 554 fbdata%ivlqcf(2,1,jo,1) = sstdata%nqc(jo) - 10 555 ELSE 556 fbdata%ivqc(jo,1) = MAX(sstdata%nqc(jo),1) 557 fbdata%ivlqc(1,jo,1) = MAX(sstdata%nqc(jo),1) 558 fbdata%ivlqcf(:,1,jo,1) = 0 559 ENDIF 560 fbdata%iobsk(1,jo,1) = 0 561 DO ja = 1, nadd 562 fbdata%padd(1,jo,1+ja,1) = & 563 & sstdata%rext(jo,padd%ipoint(ja)) 564 END DO 565 DO je = 1, next 566 fbdata%pext(1,jo,je) = & 567 & sstdata%rext(jo,pext%ipoint(je)) 568 END DO 569 570 END DO 571 572 ! Write the obfbdata structure 573 574 CALL write_obfbdata( cfname, fbdata ) 575 576 ! Output some basic statistics 577 CALL obs_wri_stats( fbdata ) 578 579 CALL dealloc_obfbdata( fbdata ) 580 581 END SUBROUTINE obs_wri_sst 582 583 SUBROUTINE obs_wri_sss 584 END SUBROUTINE obs_wri_sss 585 586 SUBROUTINE obs_wri_seaice( cprefix, seaicedata, padd, pext ) 587 !!----------------------------------------------------------------------- 588 !! 589 !! *** ROUTINE obs_wri_seaice *** 590 !! 591 !! ** Purpose : Write sea ice observation diagnostics 592 !! related 593 !! 594 !! ** Method : NetCDF 595 !! 596 !! ** Action : 597 !! 598 !! ! 07-07 (S. Ricci) Original 599 !! ! 09-01 (K. Mogensen) New feedback format. 600 !!----------------------------------------------------------------------- 601 602 !! * Modules used 603 IMPLICIT NONE 604 605 !! * Arguments 606 CHARACTER(LEN=*), INTENT(IN) :: cprefix ! Prefix for output files 607 TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of sea ice 608 TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable 609 TYPE(obswriinfo), OPTIONAL :: pext ! Extra info 610 611 !! * Local declarations 612 TYPE(obfbdata) :: fbdata 613 CHARACTER(LEN=40) :: cfname ! netCDF filename 614 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_seaice' 615 INTEGER :: jo 616 INTEGER :: ja 617 INTEGER :: je 618 INTEGER :: nadd 619 INTEGER :: next 620 621 IF ( PRESENT( padd ) ) THEN 622 nadd = padd%inum 623 ELSE 624 nadd = 0 625 ENDIF 626 627 IF ( PRESENT( pext ) ) THEN 628 next = pext%inum 629 ELSE 630 next = 0 631 ENDIF 632 633 CALL init_obfbdata( fbdata ) 634 635 CALL alloc_obfbdata( fbdata, 1, seaicedata%nsurf, 1, 1, 0, .TRUE. ) 636 637 fbdata%cname(1) = 'SEAICE' 638 fbdata%coblong(1) = 'Sea ice' 639 fbdata%cobunit(1) = 'Fraction' 640 DO je = 1, next 641 fbdata%cextname(je) = pext%cdname(je) 642 fbdata%cextlong(je) = pext%cdlong(je,1) 643 fbdata%cextunit(je) = pext%cdunit(je,1) 644 END DO 645 fbdata%caddname(1) = 'Hx' 646 fbdata%caddlong(1,1) = 'Model interpolated ICE' 647 fbdata%caddunit(1,1) = 'Fraction' 648 fbdata%cgrid(1) = 'T' 649 DO ja = 1, nadd 650 fbdata%caddname(1+ja) = padd%cdname(ja) 651 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 652 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 653 END DO 654 655 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 656 657 IF(lwp) THEN 658 WRITE(numout,*) 659 WRITE(numout,*)'obs_wri_seaice :' 660 WRITE(numout,*)'~~~~~~~~~~~~~~~~' 661 WRITE(numout,*)'Writing SEAICE feedback file : ',TRIM(cfname) 662 ENDIF 663 664 ! Transform obs_prof data structure into obfbdata structure 665 fbdata%cdjuldref = '19500101000000' 666 DO jo = 1, seaicedata%nsurf 667 fbdata%plam(jo) = seaicedata%rlam(jo) 668 fbdata%pphi(jo) = seaicedata%rphi(jo) 669 WRITE(fbdata%cdtyp(jo),'(I4)') seaicedata%ntyp(jo) 670 fbdata%ivqc(jo,:) = 0 671 fbdata%ivqcf(:,jo,:) = 0 672 IF ( seaicedata%nqc(jo) > 10 ) THEN 673 fbdata%ioqc(jo) = 4 674 fbdata%ioqcf(1,jo) = 0 675 fbdata%ioqcf(2,jo) = seaicedata%nqc(jo) - 10 676 ELSE 677 fbdata%ioqc(jo) = MAX(seaicedata%nqc(jo),1) 678 fbdata%ioqcf(:,jo) = 0 679 ENDIF 680 fbdata%ipqc(jo) = 0 681 fbdata%ipqcf(:,jo) = 0 682 fbdata%itqc(jo) = 0 683 fbdata%itqcf(:,jo) = 0 684 fbdata%cdwmo(jo) = '' 685 fbdata%kindex(jo) = seaicedata%nsfil(jo) 686 IF (ln_grid_global) THEN 687 fbdata%iobsi(jo,1) = seaicedata%mi(jo) 688 fbdata%iobsj(jo,1) = seaicedata%mj(jo) 689 ELSE 690 fbdata%iobsi(jo,1) = mig(seaicedata%mi(jo)) 691 fbdata%iobsj(jo,1) = mjg(seaicedata%mj(jo)) 692 ENDIF 693 CALL greg2jul( 0, & 694 & seaicedata%nmin(jo), & 695 & seaicedata%nhou(jo), & 696 & seaicedata%nday(jo), & 697 & seaicedata%nmon(jo), & 698 & seaicedata%nyea(jo), & 699 & fbdata%ptim(jo), & 700 & krefdate = 19500101 ) 701 fbdata%padd(1,jo,1,1) = seaicedata%rmod(jo,1) 702 fbdata%pob(1,jo,1) = seaicedata%robs(jo,1) 703 fbdata%pdep(1,jo) = 0.0 704 fbdata%idqc(1,jo) = 0 705 fbdata%idqcf(:,1,jo) = 0 706 IF ( seaicedata%nqc(jo) > 10 ) THEN 707 fbdata%ivlqc(1,jo,1) = 4 708 fbdata%ivlqcf(1,1,jo,1) = 0 709 fbdata%ivlqcf(2,1,jo,1) = seaicedata%nqc(jo) - 10 710 ELSE 711 fbdata%ivlqc(1,jo,1) = MAX(seaicedata%nqc(jo),1) 712 fbdata%ivlqcf(:,1,jo,1) = 0 713 ENDIF 714 fbdata%iobsk(1,jo,1) = 0 715 DO ja = 1, nadd 716 fbdata%padd(1,jo,1+ja,1) = & 717 & seaicedata%rext(jo,padd%ipoint(ja)) 718 END DO 719 DO je = 1, next 720 fbdata%pext(1,jo,je) = & 721 & seaicedata%rext(jo,pext%ipoint(je)) 722 END DO 723 724 END DO 725 726 ! Write the obfbdata structure 727 CALL write_obfbdata( cfname, fbdata ) 728 729 ! Output some basic statistics 730 CALL obs_wri_stats( fbdata ) 731 732 CALL dealloc_obfbdata( fbdata ) 733 734 END SUBROUTINE obs_wri_seaice 735 736 SUBROUTINE obs_wri_vel( cprefix, profdata, k2dint, padd, pext ) 737 !!----------------------------------------------------------------------- 738 !! 739 !! *** ROUTINE obs_wri_vel *** 740 !! 741 !! ** Purpose : Write current (profile) observation 742 !! related diagnostics 743 !! 744 !! ** Method : NetCDF 745 !! 746 !! ** Action : 747 !! 748 !! History : 749 !! ! 09-01 (K. Mogensen) New feedback format routine 750 !!----------------------------------------------------------------------- 751 752 !! * Modules used 753 754 !! * Arguments 755 CHARACTER(LEN=*), INTENT(IN) :: cprefix ! Prefix for output files 756 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 757 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation method 758 TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable 759 TYPE(obswriinfo), OPTIONAL :: pext ! Extra info 760 761 !! * Local declarations 762 TYPE(obfbdata) :: fbdata 763 CHARACTER(LEN=40) :: cfname 764 INTEGER :: ilevel 765 INTEGER :: jvar 766 INTEGER :: jk 767 INTEGER :: ik 768 INTEGER :: jo 769 INTEGER :: ja 770 INTEGER :: je 771 INTEGER :: nadd 772 INTEGER :: next 773 REAL(wp) :: zpres 774 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 775 & zu, & 776 & zv 777 778 IF ( PRESENT( padd ) ) THEN 779 nadd = padd%inum 780 ELSE 781 nadd = 0 782 ENDIF 783 784 IF ( PRESENT( pext ) ) THEN 785 next = pext%inum 786 ELSE 787 next = 0 788 ENDIF 789 790 CALL init_obfbdata( fbdata ) 791 792 ! Find maximum level 793 ilevel = 0 794 DO jvar = 1, 2 795 ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 796 END DO 797 CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. ) 798 799 fbdata%cname(1) = 'UVEL' 800 fbdata%cname(2) = 'VVEL' 801 fbdata%coblong(1) = 'Zonal velocity' 802 fbdata%coblong(2) = 'Meridional velocity' 803 fbdata%cobunit(1) = 'm/s' 804 fbdata%cobunit(2) = 'm/s' 805 DO je = 1, next 806 fbdata%cextname(je) = pext%cdname(je) 807 fbdata%cextlong(je) = pext%cdlong(je,1) 808 fbdata%cextunit(je) = pext%cdunit(je,1) 809 END DO 810 fbdata%caddname(1) = 'Hx' 811 fbdata%caddlong(1,1) = 'Model interpolated zonal velocity' 812 fbdata%caddlong(1,2) = 'Model interpolated meridional velocity' 813 fbdata%caddunit(1,1) = 'm/s' 814 fbdata%caddunit(1,2) = 'm/s' 815 fbdata%caddname(2) = 'HxG' 816 fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)' 817 fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)' 818 fbdata%caddunit(2,1) = 'm/s' 819 fbdata%caddunit(2,2) = 'm/s' 820 fbdata%cgrid(1) = 'U' 821 fbdata%cgrid(2) = 'V' 822 DO ja = 1, nadd 823 fbdata%caddname(2+ja) = padd%cdname(ja) 824 fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 825 fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 826 END DO 827 828 WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 829 830 IF(lwp) THEN 831 WRITE(numout,*) 832 WRITE(numout,*)'obs_wri_vel :' 833 WRITE(numout,*)'~~~~~~~~~~~~~' 834 WRITE(numout,*)'Writing velocuty feedback file : ',TRIM(cfname) 835 ENDIF 836 837 ALLOCATE( & 838 & zu(profdata%nvprot(1)), & 839 & zv(profdata%nvprot(2)) & 840 & ) 841 CALL obs_rotvel( profdata, k2dint, zu, zv ) 842 843 ! Transform obs_prof data structure into obfbdata structure 844 fbdata%cdjuldref = '19500101000000' 845 DO jo = 1, profdata%nprof 846 fbdata%plam(jo) = profdata%rlam(jo) 847 fbdata%pphi(jo) = profdata%rphi(jo) 848 WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo) 849 fbdata%ivqc(jo,:) = profdata%ivqc(jo,:) 850 fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 851 IF ( profdata%nqc(jo) > 10 ) THEN 852 fbdata%ioqc(jo) = 4 853 fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 854 fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10 855 ELSE 856 fbdata%ioqc(jo) = profdata%nqc(jo) 857 fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo) 858 ENDIF 859 fbdata%ipqc(jo) = profdata%ipqc(jo) 860 fbdata%ipqcf(:,jo) = profdata%ipqcf(:,jo) 861 fbdata%itqc(jo) = profdata%itqc(jo) 862 fbdata%itqcf(:,jo) = profdata%itqcf(:,jo) 863 fbdata%cdwmo(jo) = profdata%cwmo(jo) 864 fbdata%kindex(jo) = profdata%npfil(jo) 865 DO jvar = 1, profdata%nvar 866 IF (ln_grid_global) THEN 867 fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar) 868 fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar) 869 ELSE 870 fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar)) 871 fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar)) 872 ENDIF 873 END DO 874 CALL greg2jul( 0, & 875 & profdata%nmin(jo), & 876 & profdata%nhou(jo), & 877 & profdata%nday(jo), & 878 & profdata%nmon(jo), & 879 & profdata%nyea(jo), & 880 & fbdata%ptim(jo), & 881 & krefdate = 19500101 ) 882 ! Reform the profiles arrays for output 883 DO jvar = 1, 2 884 DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 885 ik = profdata%var(jvar)%nvlidx(jk) 886 IF ( jvar == 1 ) THEN 887 fbdata%padd(ik,jo,1,jvar) = zu(jk) 888 ELSE 889 fbdata%padd(ik,jo,1,jvar) = zv(jk) 890 ENDIF 891 fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk) 892 fbdata%pob(ik,jo,jvar) = profdata%var(jvar)%vobs(jk) 893 fbdata%pdep(ik,jo) = profdata%var(jvar)%vdep(jk) 894 fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk) 895 fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk) 896 IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN 897 fbdata%ivlqc(ik,jo,jvar) = 4 898 fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 899 fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 900 ELSE 901 fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) 902 fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk) 903 ENDIF 904 fbdata%iobsk(ik,jo,jvar) = profdata%var(jvar)%mvk(jk) 905 DO ja = 1, nadd 906 fbdata%padd(ik,jo,2+ja,jvar) = & 907 & profdata%var(jvar)%vext(jk,padd%ipoint(ja)) 908 END DO 909 DO je = 1, next 910 fbdata%pext(ik,jo,je) = & 911 & profdata%var(jvar)%vext(jk,pext%ipoint(je)) 912 END DO 913 END DO 914 END DO 915 END DO 916 917 ! Write the obfbdata structure 918 CALL write_obfbdata( cfname, fbdata ) 919 920 ! Output some basic statistics 921 CALL obs_wri_stats( fbdata ) 922 923 CALL dealloc_obfbdata( fbdata ) 924 925 DEALLOCATE( & 926 & zu, & 927 & zv & 928 & ) 929 930 END SUBROUTINE obs_wri_vel 535 END SUBROUTINE obs_wri_surf 931 536 932 537 SUBROUTINE obs_wri_stats( fbdata ) … … 952 557 INTEGER :: jk 953 558 954 ! INTEGER :: nlev955 ! INTEGER :: nlevmpp956 ! INTEGER :: nobsmpp957 559 INTEGER :: numgoodobs 958 560 INTEGER :: numgoodobsmpp
Note: See TracChangeset
for help on using the changeset viewer.