- Timestamp:
- 2015-12-04T11:56:46+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r5990 r5998 7 7 8 8 !!---------------------------------------------------------------------- 9 !! obs_pro _opt : Compute the model counterpart of temperature and10 !! salinity observations from profiles9 !! obs_prof_opt : Compute the model counterpart of profile data 10 !! obs_surf_opt : Compute the model counterpart of surface data 11 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and 12 12 !! salinity observations from profiles in generalised 13 13 !! vertical coordinates 14 !! obs_sla_opt : Compute the model counterpart of sea level anomaly15 !! observations16 !! obs_sst_opt : Compute the model counterpart of sea surface temperature17 !! observations18 !! obs_sss_opt : Compute the model counterpart of sea surface salinity19 !! observations20 !! obs_seaice_opt : Compute the model counterpart of sea ice concentration21 !! observations22 !!23 !! obs_vel_opt : Compute the model counterpart of zonal and meridional24 !! components of velocity from observations.25 14 !!---------------------------------------------------------------------- 26 15 27 !! * Modules used 16 !! * Modules used 28 17 USE par_kind, ONLY : & ! Precision variables 29 18 & wp 30 19 USE in_out_manager ! I/O manager 31 20 USE obs_inter_sup ! Interpolation support 32 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs ervationpt21 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 33 22 & obs_int_h2d, & 34 23 & obs_int_h2d_init 35 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs ervationpt24 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 36 25 & obs_int_z1d, & 37 26 & obs_int_z1d_spl … … 50 39 USE obs_grid, ONLY : & 51 40 & obs_level_search 52 41 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 42 & sbc_dcy, nday_qsr 43 53 44 IMPLICIT NONE 54 45 … … 56 47 PRIVATE 57 48 58 PUBLIC obs_pro _opt, & ! Compute the model counterpart of profile observations49 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 59 50 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations 60 ! in generalised vertical coordinates 61 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 62 & obs_sst_opt, & ! Compute the model counterpart of SST observations 63 & obs_sss_opt, & ! Compute the model counterpart of SSS observations 64 & obs_seaice_opt, & 65 & obs_vel_opt ! Compute the model counterpart of velocity profile data 66 67 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 51 & obs_surf_opt ! Compute the model counterpart of surface obs 52 53 INTEGER, PARAMETER, PUBLIC :: & 54 & imaxavtypes = 20 ! Max number of daily avgd obs types 68 55 69 56 !!---------------------------------------------------------------------- … … 77 64 CONTAINS 78 65 79 SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 80 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 81 & kdailyavtypes ) 66 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 67 & kit000, kdaystp, & 68 & pvar1, pvar2, pgdept, pmask1, pmask2, & 69 & plam1, plam2, pphi1, pphi2, & 70 & k1dint, k2dint, kdailyavtypes ) 71 82 72 !!----------------------------------------------------------------------- 83 73 !! … … 92 82 !! 93 83 !! First, a vertical profile of horizontally interpolated model 94 !! now temperatures is computed at the obs (lon, lat) point.84 !! now values is computed at the obs (lon, lat) point. 95 85 !! Several horizontal interpolation schemes are available: 96 86 !! - distance-weighted (great circle) (k2dint = 0) … … 100 90 !! - polynomial (quadrilateral grid) (k2dint = 4) 101 91 !! 102 !! Next, the vertical temperatureprofile is interpolated to the92 !! Next, the vertical profile is interpolated to the 103 93 !! data depth points. Two vertical interpolation schemes are 104 94 !! available: … … 110 100 !! routine. 111 101 !! 112 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is102 !! If the logical is switched on, the model equivalent is 113 103 !! a daily mean model temperature field. So, we first compute 114 104 !! the mean, then interpolate only at the end of the day. 115 105 !! 116 !! Note: thein situ temperature observations must be converted106 !! Note: in situ temperature observations must be converted 117 107 !! to potential temperature (the model variable) prior to 118 108 !! assimilation. 119 !!??????????????????????????????????????????????????????????????120 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???121 !!??????????????????????????????????????????????????????????????122 109 !! 123 110 !! ** Action : … … 129 116 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 130 117 !! ! 07-03 (K. Mogensen) General handling of profiles 118 !! ! 15-02 (M. Martin) Combined routine for all profile types 131 119 !!----------------------------------------------------------------------- 132 120 133 121 !! * Modules used 134 122 USE obs_profiles_def ! Definition of storage space for profile obs. … … 137 125 138 126 !! * Arguments 139 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 140 INTEGER, INTENT(IN) :: kt ! Time step 141 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 127 TYPE(obs_prof), INTENT(INOUT) :: & 128 & prodatqc ! Subset of profile data passing QC 129 INTEGER, INTENT(IN) :: kt ! Time step 130 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 142 131 INTEGER, INTENT(IN) :: kpj 143 132 INTEGER, INTENT(IN) :: kpk 144 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step145 146 INTEGER, INTENT(IN) :: k1dint 147 INTEGER, INTENT(IN) :: k2dint 148 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day133 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 134 ! (kit000-1 = restart time) 135 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 136 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 137 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 149 138 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 150 & ptn, & ! Model temperature field 151 & psn, & ! Model salinity field 152 & ptmask ! Land-sea mask 139 & pvar1, & ! Model field 1 140 & pvar2, & ! Model field 2 141 & pmask1, & ! Land-sea mask 1 142 & pmask2 ! Land-sea mask 2 143 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 144 & plam1, & ! Model longitudes for variable 1 145 & plam2, & ! Model longitudes for variable 2 146 & pphi1, & ! Model latitudes for variable 1 147 & pphi2 ! Model latitudes for variable 2 153 148 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 154 & pgdept ! Model array of depth levels149 & pgdept ! Model array of depth levels 155 150 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 156 & kdailyavtypes! Types for daily averages 151 & kdailyavtypes ! Types for daily averages 152 157 153 !! * Local declarations 158 154 INTEGER :: ji … … 168 164 INTEGER, DIMENSION(imaxavtypes) :: & 169 165 & idailyavtypes 166 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 167 & igrdi1, & 168 & igrdi2, & 169 & igrdj1, & 170 & igrdj2 170 171 REAL(KIND=wp) :: zlam 171 172 REAL(KIND=wp) :: zphi 172 173 REAL(KIND=wp) :: zdaystp 173 174 REAL(KIND=wp), DIMENSION(kpk) :: & 174 & zobsmask, & 175 & zobsmask1, & 176 & zobsmask2, & 175 177 & zobsk, & 176 178 & zobs2k 177 179 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 178 & zweig 180 & zweig1, & 181 & zweig2 179 182 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 180 & zmask, & 181 & zintt, & 182 & zints, & 183 & zinmt, & 184 & zinms 183 & zmask1, & 184 & zmask2, & 185 & zint1, & 186 & zint2, & 187 & zinm1, & 188 & zinm2 185 189 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 186 & zglam , &187 & zg phi188 INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::&189 & igrdi, &190 & igrdj190 & zglam1, & 191 & zglam2, & 192 & zgphi1, & 193 & zgphi2 194 LOGICAL :: ld_dailyav 191 195 192 196 !------------------------------------------------------------------------ 193 197 ! Local initialization 194 198 !------------------------------------------------------------------------ 195 ! ...Record and data counters199 ! Record and data counters 196 200 inrc = kt - kit000 + 2 197 201 ipro = prodatqc%npstp(inrc) 198 202 199 203 ! Daily average types 204 ld_dailyav = .FALSE. 200 205 IF ( PRESENT(kdailyavtypes) ) THEN 201 206 idailyavtypes(:) = kdailyavtypes(:) 207 IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 202 208 ELSE 203 209 idailyavtypes(:) = -1 204 210 ENDIF 205 211 206 ! Initialize daily mean for first timestep 212 ! Daily means are calculated for values over timesteps: 213 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 207 214 idayend = MOD( kt - kit000 + 1, kdaystp ) 208 215 209 ! Added kt == 0 test to catch restart case 210 IF ( idayend == 1 .OR. kt == 0) THEN 211 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 216 IF ( ld_dailyav ) THEN 217 218 ! Initialize daily mean for first timestep of the day 219 IF ( idayend == 1 .OR. kt == 0 ) THEN 220 DO jk = 1, jpk 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 prodatqc%vdmean(ji,jj,jk,1) = 0.0 224 prodatqc%vdmean(ji,jj,jk,2) = 0.0 225 END DO 226 END DO 227 END DO 228 ENDIF 229 212 230 DO jk = 1, jpk 213 231 DO jj = 1, jpj 214 232 DO ji = 1, jpi 215 prodatqc%vdmean(ji,jj,jk,1) = 0.0 216 prodatqc%vdmean(ji,jj,jk,2) = 0.0 233 ! Increment field 1 for computing daily mean 234 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 235 & + pvar1(ji,jj,jk) 236 ! Increment field 2 for computing daily mean 237 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 238 & + pvar2(ji,jj,jk) 217 239 END DO 218 240 END DO 219 241 END DO 220 ENDIF 221 222 DO jk = 1, jpk 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 ! Increment the temperature field for computing daily mean 226 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 227 & + ptn(ji,jj,jk) 228 ! Increment the salinity field for computing daily mean 229 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 230 & + psn(ji,jj,jk) 231 END DO 232 END DO 233 END DO 234 235 ! Compute the daily mean at the end of day 236 zdaystp = 1.0 / REAL( kdaystp ) 237 IF ( idayend == 0 ) THEN 238 DO jk = 1, jpk 239 DO jj = 1, jpj 240 DO ji = 1, jpi 241 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 242 & * zdaystp 243 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 244 & * zdaystp 242 243 ! Compute the daily mean at the end of day 244 zdaystp = 1.0 / REAL( kdaystp ) 245 IF ( idayend == 0 ) THEN 246 IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 247 CALL FLUSH(numout) 248 DO jk = 1, jpk 249 DO jj = 1, jpj 250 DO ji = 1, jpi 251 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 252 & * zdaystp 253 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 254 & * zdaystp 255 END DO 245 256 END DO 246 257 END DO 247 END DO 258 ENDIF 259 248 260 ENDIF 249 261 250 262 ! Get the data for interpolation 251 263 ALLOCATE( & 252 & igrdi(2,2,ipro), & 253 & igrdj(2,2,ipro), & 254 & zglam(2,2,ipro), & 255 & zgphi(2,2,ipro), & 256 & zmask(2,2,kpk,ipro), & 257 & zintt(2,2,kpk,ipro), & 258 & zints(2,2,kpk,ipro) & 264 & igrdi1(2,2,ipro), & 265 & igrdi2(2,2,ipro), & 266 & igrdj1(2,2,ipro), & 267 & igrdj2(2,2,ipro), & 268 & zglam1(2,2,ipro), & 269 & zglam2(2,2,ipro), & 270 & zgphi1(2,2,ipro), & 271 & zgphi2(2,2,ipro), & 272 & zmask1(2,2,kpk,ipro), & 273 & zmask2(2,2,kpk,ipro), & 274 & zint1(2,2,kpk,ipro), & 275 & zint2(2,2,kpk,ipro) & 259 276 & ) 260 277 261 278 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 262 279 iobs = jobs - prodatqc%nprofup 263 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 264 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 265 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 266 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 267 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 268 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 269 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 270 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 280 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 281 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 282 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 283 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 284 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 285 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 286 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 287 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 288 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 289 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 290 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 291 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 292 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 293 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 294 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 295 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 271 296 END DO 272 297 273 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 274 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 275 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 276 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 277 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 298 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 299 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 300 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 301 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 302 303 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 304 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 305 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 306 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 278 307 279 308 ! At the end of the day also get interpolated means 280 IF ( idayend == 0 ) THEN309 IF ( ld_dailyav .AND. idayend == 0 ) THEN 281 310 282 311 ALLOCATE( & 283 & zinm t(2,2,kpk,ipro), &284 & zinm s(2,2,kpk,ipro) &312 & zinm1(2,2,kpk,ipro), & 313 & zinm2(2,2,kpk,ipro) & 285 314 & ) 286 315 287 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &288 & prodatqc%vdmean(:,:,:,1), zinm t)289 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &290 & prodatqc%vdmean(:,:,:,2), zinm s)316 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 317 & prodatqc%vdmean(:,:,:,1), zinm1 ) 318 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 319 & prodatqc%vdmean(:,:,:,2), zinm2 ) 291 320 292 321 ENDIF … … 297 326 298 327 IF ( kt /= prodatqc%mstp(jobs) ) THEN 299 328 300 329 IF(lwp) THEN 301 330 WRITE(numout,*) … … 312 341 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 313 342 ENDIF 314 343 315 344 zlam = prodatqc%rlam(jobs) 316 345 zphi = prodatqc%rphi(jobs) 317 346 318 347 ! Horizontal weights and vertical mask 319 348 320 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 321 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 349 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 322 350 323 351 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 324 & zglam (:,:,iobs), zgphi(:,:,iobs), &325 & zmask (:,:,:,iobs), zweig, zobsmask)352 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 353 & zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 326 354 327 355 ENDIF 328 356 357 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 358 359 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 360 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 361 & zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 362 363 ENDIF 364 329 365 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 330 366 331 367 zobsk(:) = obfillflt 332 368 333 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN369 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 334 370 335 371 IF ( idayend == 0 ) THEN 336 337 ! Daily averaged moored buoy (MRB) data 338 372 ! Daily averaged data 339 373 CALL obs_int_h2d( kpk, kpk, & 340 & zweig, zinmt(:,:,:,iobs), zobsk ) 341 342 343 ELSE 344 345 CALL ctl_stop( ' A nonzero' // & 346 & ' number of profile T BUOY data should' // & 347 & ' only occur at the end of a given day' ) 374 & zweig1, zinm1(:,:,:,iobs), zobsk ) 348 375 349 376 ENDIF 350 377 351 378 ELSE 352 379 353 380 ! Point data 354 355 381 CALL obs_int_h2d( kpk, kpk, & 356 & zweig , zintt(:,:,:,iobs), zobsk )382 & zweig1, zint1(:,:,:,iobs), zobsk ) 357 383 358 384 ENDIF … … 362 388 ! polynomial at obs points 363 389 !------------------------------------------------------------- 364 390 365 391 IF ( k1dint == 1 ) THEN 366 392 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 367 & pgdept, zobsmask )393 & pgdept, zobsmask1 ) 368 394 ENDIF 369 395 370 396 !----------------------------------------------------------------- 371 397 ! Vertical interpolation to the observation point … … 379 405 & zobsk, zobs2k, & 380 406 & prodatqc%var(1)%vmod(ista:iend), & 381 & pgdept, zobsmask )407 & pgdept, zobsmask1 ) 382 408 383 409 ENDIF … … 391 417 IF ( idayend == 0 ) THEN 392 418 393 ! Daily averaged moored buoy (MRB) data 394 419 ! Daily averaged data 395 420 CALL obs_int_h2d( kpk, kpk, & 396 & zweig, zinms(:,:,:,iobs), zobsk ) 397 398 ELSE 399 400 CALL ctl_stop( ' A nonzero' // & 401 & ' number of profile S BUOY data should' // & 402 & ' only occur at the end of a given day' ) 421 & zweig2, zinm2(:,:,:,iobs), zobsk ) 403 422 404 423 ENDIF 405 424 406 425 ELSE 407 426 408 427 ! Point data 409 410 428 CALL obs_int_h2d( kpk, kpk, & 411 & zweig , zints(:,:,:,iobs), zobsk )429 & zweig2, zint2(:,:,:,iobs), zobsk ) 412 430 413 431 ENDIF … … 418 436 ! polynomial at obs points 419 437 !------------------------------------------------------------- 420 438 421 439 IF ( k1dint == 1 ) THEN 422 440 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 423 & pgdept, zobsmask )441 & pgdept, zobsmask2 ) 424 442 ENDIF 425 443 426 444 !---------------------------------------------------------------- 427 445 ! Vertical interpolation to the observation point … … 435 453 & zobsk, zobs2k, & 436 454 & prodatqc%var(2)%vmod(ista:iend),& 437 & pgdept, zobsmask )455 & pgdept, zobsmask2 ) 438 456 439 457 ENDIF 440 458 441 459 END DO 442 460 443 461 ! Deallocate the data for interpolation 444 462 DEALLOCATE( & 445 & igrdi, & 446 & igrdj, & 447 & zglam, & 448 & zgphi, & 449 & zmask, & 450 & zintt, & 451 & zints & 463 & igrdi1, & 464 & igrdi2, & 465 & igrdj1, & 466 & igrdj2, & 467 & zglam1, & 468 & zglam2, & 469 & zgphi1, & 470 & zgphi2, & 471 & zmask1, & 472 & zmask2, & 473 & zint1, & 474 & zint2 & 452 475 & ) 476 453 477 ! At the end of the day also get interpolated means 454 IF ( idayend == 0 ) THEN478 IF ( ld_dailyav .AND. idayend == 0 ) THEN 455 479 DEALLOCATE( & 456 & zinm t, &457 & zinm s&480 & zinm1, & 481 & zinm2 & 458 482 & ) 459 483 ENDIF 460 484 461 485 prodatqc%nprofup = prodatqc%nprofup + ipro 462 463 END SUBROUTINE obs_pro _opt486 487 END SUBROUTINE obs_prof_opt 464 488 465 489 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & … … 1054 1078 END SUBROUTINE obs_pro_sco_opt 1055 1079 1056 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 1057 & psshn, psshmask, k2dint ) 1080 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 1081 & kit000, kdaystp, psurf, psurfmask, & 1082 & k2dint, ldnightav ) 1083 1058 1084 !!----------------------------------------------------------------------- 1059 1085 !! 1060 !! *** ROUTINE obs_s la_opt ***1061 !! 1062 !! ** Purpose : Compute the model counterpart of s ea level anomaly1086 !! *** ROUTINE obs_surf_opt *** 1087 !! 1088 !! ** Purpose : Compute the model counterpart of surface 1063 1089 !! data by interpolating from the model grid to the 1064 1090 !! observation point. … … 1067 1093 !! the model values at the corners of the surrounding grid box. 1068 1094 !! 1069 !! The n ow model SSHis first computed at the obs (lon, lat) point.1095 !! The new model value is first computed at the obs (lon, lat) point. 1070 1096 !! 1071 1097 !! Several horizontal interpolation schemes are available: … … 1075 1101 !! - bilinear (quadrilateral grid) (k2dint = 3) 1076 1102 !! - polynomial (quadrilateral grid) (k2dint = 4) 1077 !! 1078 !! The sea level anomaly at the observation points is then computed 1079 !! by removing a mean dynamic topography (defined at the obs. point). 1103 !! 1080 1104 !! 1081 1105 !! ** Action : … … 1083 1107 !! History : 1084 1108 !! ! 07-03 (A. Weaver) 1109 !! ! 15-02 (M. Martin) Combined routine for surface types 1085 1110 !!----------------------------------------------------------------------- 1086 1111 1087 1112 !! * Modules used 1088 1113 USE obs_surf_def ! Definition of storage space for surface observations … … 1091 1116 1092 1117 !! * Arguments 1093 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening 1094 INTEGER, INTENT(IN) :: kt ! Time step 1095 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1118 TYPE(obs_surf), INTENT(INOUT) :: & 1119 & surfdataqc ! Subset of surface data passing QC 1120 INTEGER, INTENT(IN) :: kt ! Time step 1121 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1096 1122 INTEGER, INTENT(IN) :: kpj 1097 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1098 ! (kit000-1 = restart time) 1099 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1100 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 1101 & psshn, & ! Model SSH field 1102 & psshmask ! Land-sea mask 1103 1123 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1124 ! (kit000-1 = restart time) 1125 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 1126 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1127 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 1128 & psurf, & ! Model surface field 1129 & psurfmask ! Land-sea mask 1130 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 1131 1104 1132 !! * Local declarations 1105 1133 INTEGER :: ji … … 1107 1135 INTEGER :: jobs 1108 1136 INTEGER :: inrc 1109 INTEGER :: is la1137 INTEGER :: isurf 1110 1138 INTEGER :: iobs 1111 REAL(KIND=wp) :: zlam 1112 REAL(KIND=wp) :: zphi 1113 REAL(KIND=wp) :: zext(1), zobsmask(1) 1114 REAL(kind=wp), DIMENSION(2,2,1) :: & 1115 & zweig 1116 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1117 & zmask, & 1118 & zsshl, & 1119 & zglam, & 1120 & zgphi 1139 INTEGER :: idayend 1121 1140 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1122 1141 & igrdi, & 1123 1142 & igrdj 1143 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1144 & icount_night, & 1145 & imask_night 1146 REAL(wp) :: zlam 1147 REAL(wp) :: zphi 1148 REAL(wp), DIMENSION(1) :: zext, zobsmask 1149 REAL(wp) :: zdaystp 1150 REAL(wp), DIMENSION(2,2,1) :: & 1151 & zweig 1152 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1153 & zmask, & 1154 & zsurf, & 1155 & zsurfm, & 1156 & zglam, & 1157 & zgphi 1158 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1159 & zintmp, & 1160 & zouttmp, & 1161 & zmeanday ! to compute model sst in region of 24h daylight (pole) 1124 1162 1125 1163 !------------------------------------------------------------------------ 1126 1164 ! Local initialization 1127 1165 !------------------------------------------------------------------------ 1128 ! ...Record and data counters1166 ! Record and data counters 1129 1167 inrc = kt - kit000 + 2 1130 isla = sladatqc%nsstp(inrc) 1168 isurf = surfdataqc%nsstp(inrc) 1169 1170 IF ( ldnightav ) THEN 1171 1172 ! Initialize array for night mean 1173 IF ( kt == 0 ) THEN 1174 ALLOCATE ( icount_night(kpi,kpj) ) 1175 ALLOCATE ( imask_night(kpi,kpj) ) 1176 ALLOCATE ( zintmp(kpi,kpj) ) 1177 ALLOCATE ( zouttmp(kpi,kpj) ) 1178 ALLOCATE ( zmeanday(kpi,kpj) ) 1179 nday_qsr = -1 ! initialisation flag for nbc_dcy 1180 ENDIF 1181 1182 ! Night-time means are calculated for night-time values over timesteps: 1183 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 1184 idayend = MOD( kt - kit000 + 1, kdaystp ) 1185 1186 ! Initialize night-time mean for first timestep of the day 1187 IF ( idayend == 1 .OR. kt == 0 ) THEN 1188 DO jj = 1, jpj 1189 DO ji = 1, jpi 1190 surfdataqc%vdmean(ji,jj) = 0.0 1191 zmeanday(ji,jj) = 0.0 1192 icount_night(ji,jj) = 0 1193 END DO 1194 END DO 1195 ENDIF 1196 1197 zintmp(:,:) = 0.0 1198 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 1199 imask_night(:,:) = INT( zouttmp(:,:) ) 1200 1201 DO jj = 1, jpj 1202 DO ji = 1, jpi 1203 ! Increment the temperature field for computing night mean and counter 1204 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 1205 & + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 1206 zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj) 1207 icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 1208 END DO 1209 END DO 1210 1211 ! Compute the night-time mean at the end of the day 1212 zdaystp = 1.0 / REAL( kdaystp ) 1213 IF ( idayend == 0 ) THEN 1214 IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 1215 DO jj = 1, jpj 1216 DO ji = 1, jpi 1217 ! Test if "no night" point 1218 IF ( icount_night(ji,jj) > 0 ) THEN 1219 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 1220 & / REAL( icount_night(ji,jj) ) 1221 ELSE 1222 !At locations where there is no night (e.g. poles), 1223 ! calculate daily mean instead of night-time mean. 1224 surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 1225 ENDIF 1226 END DO 1227 END DO 1228 ENDIF 1229 1230 ENDIF 1131 1231 1132 1232 ! Get the data for interpolation 1133 1233 1134 1234 ALLOCATE( & 1135 & igrdi(2,2,is la), &1136 & igrdj(2,2,is la), &1137 & zglam(2,2,is la), &1138 & zgphi(2,2,is la), &1139 & zmask(2,2,is la), &1140 & zs shl(2,2,isla) &1235 & igrdi(2,2,isurf), & 1236 & igrdj(2,2,isurf), & 1237 & zglam(2,2,isurf), & 1238 & zgphi(2,2,isurf), & 1239 & zmask(2,2,isurf), & 1240 & zsurf(2,2,isurf) & 1141 1241 & ) 1142 1143 DO jobs = s ladatqc%nsurfup + 1, sladatqc%nsurfup + isla1144 iobs = jobs - s ladatqc%nsurfup1145 igrdi(1,1,iobs) = s ladatqc%mi(jobs)-11146 igrdj(1,1,iobs) = s ladatqc%mj(jobs)-11147 igrdi(1,2,iobs) = s ladatqc%mi(jobs)-11148 igrdj(1,2,iobs) = s ladatqc%mj(jobs)1149 igrdi(2,1,iobs) = s ladatqc%mi(jobs)1150 igrdj(2,1,iobs) = s ladatqc%mj(jobs)-11151 igrdi(2,2,iobs) = s ladatqc%mi(jobs)1152 igrdj(2,2,iobs) = s ladatqc%mj(jobs)1242 1243 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1244 iobs = jobs - surfdataqc%nsurfup 1245 igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 1246 igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 1247 igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 1248 igrdj(1,2,iobs) = surfdataqc%mj(jobs) 1249 igrdi(2,1,iobs) = surfdataqc%mi(jobs) 1250 igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 1251 igrdi(2,2,iobs) = surfdataqc%mi(jobs) 1252 igrdj(2,2,iobs) = surfdataqc%mj(jobs) 1153 1253 END DO 1154 1254 1155 CALL obs_int_comm_2d( 2, 2, is la, &1255 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 1156 1256 & igrdi, igrdj, glamt, zglam ) 1157 CALL obs_int_comm_2d( 2, 2, is la, &1257 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 1158 1258 & igrdi, igrdj, gphit, zgphi ) 1159 CALL obs_int_comm_2d( 2, 2, isla, & 1160 & igrdi, igrdj, psshmask, zmask ) 1161 CALL obs_int_comm_2d( 2, 2, isla, & 1162 & igrdi, igrdj, psshn, zsshl ) 1259 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 1260 & igrdi, igrdj, psurfmask, zmask ) 1261 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 1262 & igrdi, igrdj, psurf, zsurf ) 1263 1264 ! At the end of the day get interpolated means 1265 IF ( idayend == 0 .AND. ldnightav ) THEN 1266 1267 ALLOCATE( & 1268 & zsurfm(2,2,isurf) & 1269 & ) 1270 1271 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 1272 & surfdataqc%vdmean(:,:), zsurfm ) 1273 1274 ENDIF 1163 1275 1164 1276 ! Loop over observations 1165 1166 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 1167 1168 iobs = jobs - sladatqc%nsurfup 1169 1170 IF ( kt /= sladatqc%mstp(jobs) ) THEN 1171 1277 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1278 1279 iobs = jobs - surfdataqc%nsurfup 1280 1281 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 1282 1172 1283 IF(lwp) THEN 1173 1284 WRITE(numout,*) … … 1179 1290 WRITE(numout,*) ' Record = ', jobs, & 1180 1291 & ' kt = ', kt, & 1181 & ' mstp = ', s ladatqc%mstp(jobs), &1182 & ' ntyp = ', s ladatqc%ntyp(jobs)1292 & ' mstp = ', surfdataqc%mstp(jobs), & 1293 & ' ntyp = ', surfdataqc%ntyp(jobs) 1183 1294 ENDIF 1184 CALL ctl_stop( 'obs_s la_opt', 'Inconsistent time' )1185 1295 CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 1296 1186 1297 ENDIF 1187 1188 zlam = s ladatqc%rlam(jobs)1189 zphi = s ladatqc%rphi(jobs)1190 1191 ! Get weights to interpolate the model SSHto the observation point1298 1299 zlam = surfdataqc%rlam(jobs) 1300 zphi = surfdataqc%rphi(jobs) 1301 1302 ! Get weights to interpolate the model value to the observation point 1192 1303 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 1193 1304 & zglam(:,:,iobs), zgphi(:,:,iobs), & 1194 1305 & zmask(:,:,iobs), zweig, zobsmask ) 1195 1196 1197 ! Interpolate the model SSH to the observation point 1198 CALL obs_int_h2d( 1, 1, & 1199 & zweig, zsshl(:,:,iobs), zext ) 1200 1201 sladatqc%rext(jobs,1) = zext(1) 1202 ! ... Remove the MDT at the observation point 1203 sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 1306 1307 ! Interpolate the model field to the observation point 1308 IF ( ldnightav .AND. idayend == 0 ) THEN 1309 ! Night-time averaged data 1310 CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 1311 ELSE 1312 CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext ) 1313 ENDIF 1314 1315 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 1316 ! ... Remove the MDT from the SSH at the observation point to get the SLA 1317 surfdataqc%rext(jobs,1) = zext(1) 1318 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 1319 ELSE 1320 surfdataqc%rmod(jobs,1) = zext(1) 1321 ENDIF 1204 1322 1205 1323 END DO … … 1212 1330 & zgphi, & 1213 1331 & zmask, & 1214 & zs shl&1332 & zsurf & 1215 1333 & ) 1216 1334 1217 sladatqc%nsurfup = sladatqc%nsurfup + isla 1218 1219 END SUBROUTINE obs_sla_opt 1220 1221 SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 1222 & psstn, psstmask, k2dint, ld_nightav ) 1223 !!----------------------------------------------------------------------- 1224 !! 1225 !! *** ROUTINE obs_sst_opt *** 1226 !! 1227 !! ** Purpose : Compute the model counterpart of surface temperature 1228 !! data by interpolating from the model grid to the 1229 !! observation point. 1230 !! 1231 !! ** Method : Linearly interpolate to each observation point using 1232 !! the model values at the corners of the surrounding grid box. 1233 !! 1234 !! The now model SST is first computed at the obs (lon, lat) point. 1235 !! 1236 !! Several horizontal interpolation schemes are available: 1237 !! - distance-weighted (great circle) (k2dint = 0) 1238 !! - distance-weighted (small angle) (k2dint = 1) 1239 !! - bilinear (geographical grid) (k2dint = 2) 1240 !! - bilinear (quadrilateral grid) (k2dint = 3) 1241 !! - polynomial (quadrilateral grid) (k2dint = 4) 1242 !! 1243 !! 1244 !! ** Action : 1245 !! 1246 !! History : 1247 !! ! 07-07 (S. Ricci ) : Original 1248 !! 1249 !!----------------------------------------------------------------------- 1250 1251 !! * Modules used 1252 USE obs_surf_def ! Definition of storage space for surface observations 1253 USE sbcdcy 1254 1255 IMPLICIT NONE 1256 1257 !! * Arguments 1258 TYPE(obs_surf), INTENT(INOUT) :: & 1259 & sstdatqc ! Subset of surface data not failing screening 1260 INTEGER, INTENT(IN) :: kt ! Time step 1261 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1262 INTEGER, INTENT(IN) :: kpj 1263 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1264 ! (kit000-1 = restart time) 1265 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1266 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 1267 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 1268 & psstn, & ! Model SST field 1269 & psstmask ! Land-sea mask 1270 1271 !! * Local declarations 1272 INTEGER :: ji 1273 INTEGER :: jj 1274 INTEGER :: jobs 1275 INTEGER :: inrc 1276 INTEGER :: isst 1277 INTEGER :: iobs 1278 INTEGER :: idayend 1279 REAL(KIND=wp) :: zlam 1280 REAL(KIND=wp) :: zphi 1281 REAL(KIND=wp) :: zext(1), zobsmask(1) 1282 REAL(KIND=wp) :: zdaystp 1283 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1284 & icount_sstnight, & 1285 & imask_night 1286 REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1287 & zintmp, & 1288 & zouttmp, & 1289 & zmeanday ! to compute model sst in region of 24h daylight (pole) 1290 REAL(kind=wp), DIMENSION(2,2,1) :: & 1291 & zweig 1292 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1293 & zmask, & 1294 & zsstl, & 1295 & zsstm, & 1296 & zglam, & 1297 & zgphi 1298 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1299 & igrdi, & 1300 & igrdj 1301 LOGICAL, INTENT(IN) :: ld_nightav 1302 1303 !----------------------------------------------------------------------- 1304 ! Local initialization 1305 !----------------------------------------------------------------------- 1306 ! ... Record and data counters 1307 inrc = kt - kit000 + 2 1308 isst = sstdatqc%nsstp(inrc) 1309 1310 IF ( ld_nightav ) THEN 1311 1312 ! Initialize array for night mean 1313 1314 IF ( kt .EQ. 0 ) THEN 1315 ALLOCATE ( icount_sstnight(kpi,kpj) ) 1316 ALLOCATE ( imask_night(kpi,kpj) ) 1317 ALLOCATE ( zintmp(kpi,kpj) ) 1318 ALLOCATE ( zouttmp(kpi,kpj) ) 1319 ALLOCATE ( zmeanday(kpi,kpj) ) 1320 nday_qsr = -1 ! initialisation flag for nbc_dcy 1321 ENDIF 1322 1323 ! Initialize daily mean for first timestep 1324 idayend = MOD( kt - kit000 + 1, kdaystp ) 1325 1326 ! Added kt == 0 test to catch restart case 1327 IF ( idayend == 1 .OR. kt == 0) THEN 1328 IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 1329 DO jj = 1, jpj 1330 DO ji = 1, jpi 1331 sstdatqc%vdmean(ji,jj) = 0.0 1332 zmeanday(ji,jj) = 0.0 1333 icount_sstnight(ji,jj) = 0 1334 END DO 1335 END DO 1336 ENDIF 1337 1338 zintmp(:,:) = 0.0 1339 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 1340 imask_night(:,:) = INT( zouttmp(:,:) ) 1341 1342 DO jj = 1, jpj 1343 DO ji = 1, jpi 1344 ! Increment the temperature field for computing night mean and counter 1345 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 1346 & + psstn(ji,jj)*imask_night(ji,jj) 1347 zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj) 1348 icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 1349 END DO 1350 END DO 1351 1352 ! Compute the daily mean at the end of day 1353 1354 zdaystp = 1.0 / REAL( kdaystp ) 1355 1356 IF ( idayend == 0 ) THEN 1357 DO jj = 1, jpj 1358 DO ji = 1, jpi 1359 ! Test if "no night" point 1360 IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 1361 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 1362 & / icount_sstnight(ji,jj) 1363 ELSE 1364 sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 1365 ENDIF 1366 END DO 1367 END DO 1368 ENDIF 1369 1370 ENDIF 1371 1372 ! Get the data for interpolation 1373 1374 ALLOCATE( & 1375 & igrdi(2,2,isst), & 1376 & igrdj(2,2,isst), & 1377 & zglam(2,2,isst), & 1378 & zgphi(2,2,isst), & 1379 & zmask(2,2,isst), & 1380 & zsstl(2,2,isst) & 1381 & ) 1382 1383 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 1384 iobs = jobs - sstdatqc%nsurfup 1385 igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 1386 igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 1387 igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 1388 igrdj(1,2,iobs) = sstdatqc%mj(jobs) 1389 igrdi(2,1,iobs) = sstdatqc%mi(jobs) 1390 igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 1391 igrdi(2,2,iobs) = sstdatqc%mi(jobs) 1392 igrdj(2,2,iobs) = sstdatqc%mj(jobs) 1393 END DO 1394 1395 CALL obs_int_comm_2d( 2, 2, isst, & 1396 & igrdi, igrdj, glamt, zglam ) 1397 CALL obs_int_comm_2d( 2, 2, isst, & 1398 & igrdi, igrdj, gphit, zgphi ) 1399 CALL obs_int_comm_2d( 2, 2, isst, & 1400 & igrdi, igrdj, psstmask, zmask ) 1401 CALL obs_int_comm_2d( 2, 2, isst, & 1402 & igrdi, igrdj, psstn, zsstl ) 1403 1404 ! At the end of the day get interpolated means 1405 IF ( idayend == 0 .AND. ld_nightav ) THEN 1406 1407 ALLOCATE( & 1408 & zsstm(2,2,isst) & 1409 & ) 1410 1411 CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 1412 & sstdatqc%vdmean(:,:), zsstm ) 1413 1414 ENDIF 1415 1416 ! Loop over observations 1417 1418 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 1419 1420 iobs = jobs - sstdatqc%nsurfup 1421 1422 IF ( kt /= sstdatqc%mstp(jobs) ) THEN 1423 1424 IF(lwp) THEN 1425 WRITE(numout,*) 1426 WRITE(numout,*) ' E R R O R : Observation', & 1427 & ' time step is not consistent with the', & 1428 & ' model time step' 1429 WRITE(numout,*) ' =========' 1430 WRITE(numout,*) 1431 WRITE(numout,*) ' Record = ', jobs, & 1432 & ' kt = ', kt, & 1433 & ' mstp = ', sstdatqc%mstp(jobs), & 1434 & ' ntyp = ', sstdatqc%ntyp(jobs) 1435 ENDIF 1436 CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 1437 1438 ENDIF 1439 1440 zlam = sstdatqc%rlam(jobs) 1441 zphi = sstdatqc%rphi(jobs) 1442 1443 ! Get weights to interpolate the model SST to the observation point 1444 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 1445 & zglam(:,:,iobs), zgphi(:,:,iobs), & 1446 & zmask(:,:,iobs), zweig, zobsmask ) 1447 1448 ! Interpolate the model SST to the observation point 1449 1450 IF ( ld_nightav ) THEN 1451 1452 IF ( idayend == 0 ) THEN 1453 ! Daily averaged/diurnal cycle of SST data 1454 CALL obs_int_h2d( 1, 1, & 1455 & zweig, zsstm(:,:,iobs), zext ) 1456 ELSE 1457 CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // & 1458 & ' number of night SST data should' // & 1459 & ' only occur at the end of a given day' ) 1460 ENDIF 1461 1462 ELSE 1463 1464 CALL obs_int_h2d( 1, 1, & 1465 & zweig, zsstl(:,:,iobs), zext ) 1466 1467 ENDIF 1468 sstdatqc%rmod(jobs,1) = zext(1) 1469 1470 END DO 1471 1472 ! Deallocate the data for interpolation 1473 DEALLOCATE( & 1474 & igrdi, & 1475 & igrdj, & 1476 & zglam, & 1477 & zgphi, & 1478 & zmask, & 1479 & zsstl & 1480 & ) 1481 1482 ! At the end of the day also get interpolated means 1483 IF ( idayend == 0 .AND. ld_nightav ) THEN 1335 ! At the end of the day also deallocate night-time mean array 1336 IF ( idayend == 0 .AND. ldnightav ) THEN 1484 1337 DEALLOCATE( & 1485 & zs stm &1338 & zsurfm & 1486 1339 & ) 1487 1340 ENDIF 1488 1489 sstdatqc%nsurfup = sstdatqc%nsurfup + isst 1490 1491 END SUBROUTINE obs_sst_opt 1492 1493 SUBROUTINE obs_sss_opt 1494 !!----------------------------------------------------------------------- 1495 !! 1496 !! *** ROUTINE obs_sss_opt *** 1497 !! 1498 !! ** Purpose : Compute the model counterpart of sea surface salinity 1499 !! data by interpolating from the model grid to the 1500 !! observation point. 1501 !! 1502 !! ** Method : 1503 !! 1504 !! ** Action : 1505 !! 1506 !! History : 1507 !! ! ??-?? 1508 !!----------------------------------------------------------------------- 1509 1510 IMPLICIT NONE 1511 1512 END SUBROUTINE obs_sss_opt 1513 1514 SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, & 1515 & pseaicen, pseaicemask, k2dint ) 1516 1517 !!----------------------------------------------------------------------- 1518 !! 1519 !! *** ROUTINE obs_seaice_opt *** 1520 !! 1521 !! ** Purpose : Compute the model counterpart of surface temperature 1522 !! data by interpolating from the model grid to the 1523 !! observation point. 1524 !! 1525 !! ** Method : Linearly interpolate to each observation point using 1526 !! the model values at the corners of the surrounding grid box. 1527 !! 1528 !! The now model sea ice is first computed at the obs (lon, lat) point. 1529 !! 1530 !! Several horizontal interpolation schemes are available: 1531 !! - distance-weighted (great circle) (k2dint = 0) 1532 !! - distance-weighted (small angle) (k2dint = 1) 1533 !! - bilinear (geographical grid) (k2dint = 2) 1534 !! - bilinear (quadrilateral grid) (k2dint = 3) 1535 !! - polynomial (quadrilateral grid) (k2dint = 4) 1536 !! 1537 !! 1538 !! ** Action : 1539 !! 1540 !! History : 1541 !! ! 07-07 (S. Ricci ) : Original 1542 !! 1543 !!----------------------------------------------------------------------- 1544 1545 !! * Modules used 1546 USE obs_surf_def ! Definition of storage space for surface observations 1547 1548 IMPLICIT NONE 1549 1550 !! * Arguments 1551 TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening 1552 INTEGER, INTENT(IN) :: kt ! Time step 1553 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1554 INTEGER, INTENT(IN) :: kpj 1555 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1556 ! (kit000-1 = restart time) 1557 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1558 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 1559 & pseaicen, & ! Model sea ice field 1560 & pseaicemask ! Land-sea mask 1561 1562 !! * Local declarations 1563 INTEGER :: ji 1564 INTEGER :: jj 1565 INTEGER :: jobs 1566 INTEGER :: inrc 1567 INTEGER :: iseaice 1568 INTEGER :: iobs 1569 1570 REAL(KIND=wp) :: zlam 1571 REAL(KIND=wp) :: zphi 1572 REAL(KIND=wp) :: zext(1), zobsmask(1) 1573 REAL(kind=wp), DIMENSION(2,2,1) :: & 1574 & zweig 1575 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1576 & zmask, & 1577 & zseaicel, & 1578 & zglam, & 1579 & zgphi 1580 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1581 & igrdi, & 1582 & igrdj 1583 1584 !------------------------------------------------------------------------ 1585 ! Local initialization 1586 !------------------------------------------------------------------------ 1587 ! ... Record and data counters 1588 inrc = kt - kit000 + 2 1589 iseaice = seaicedatqc%nsstp(inrc) 1590 1591 ! Get the data for interpolation 1592 1593 ALLOCATE( & 1594 & igrdi(2,2,iseaice), & 1595 & igrdj(2,2,iseaice), & 1596 & zglam(2,2,iseaice), & 1597 & zgphi(2,2,iseaice), & 1598 & zmask(2,2,iseaice), & 1599 & zseaicel(2,2,iseaice) & 1600 & ) 1601 1602 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 1603 iobs = jobs - seaicedatqc%nsurfup 1604 igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1 1605 igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1 1606 igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1 1607 igrdj(1,2,iobs) = seaicedatqc%mj(jobs) 1608 igrdi(2,1,iobs) = seaicedatqc%mi(jobs) 1609 igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1 1610 igrdi(2,2,iobs) = seaicedatqc%mi(jobs) 1611 igrdj(2,2,iobs) = seaicedatqc%mj(jobs) 1612 END DO 1613 1614 CALL obs_int_comm_2d( 2, 2, iseaice, & 1615 & igrdi, igrdj, glamt, zglam ) 1616 CALL obs_int_comm_2d( 2, 2, iseaice, & 1617 & igrdi, igrdj, gphit, zgphi ) 1618 CALL obs_int_comm_2d( 2, 2, iseaice, & 1619 & igrdi, igrdj, pseaicemask, zmask ) 1620 CALL obs_int_comm_2d( 2, 2, iseaice, & 1621 & igrdi, igrdj, pseaicen, zseaicel ) 1622 1623 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 1624 1625 iobs = jobs - seaicedatqc%nsurfup 1626 1627 IF ( kt /= seaicedatqc%mstp(jobs) ) THEN 1628 1629 IF(lwp) THEN 1630 WRITE(numout,*) 1631 WRITE(numout,*) ' E R R O R : Observation', & 1632 & ' time step is not consistent with the', & 1633 & ' model time step' 1634 WRITE(numout,*) ' =========' 1635 WRITE(numout,*) 1636 WRITE(numout,*) ' Record = ', jobs, & 1637 & ' kt = ', kt, & 1638 & ' mstp = ', seaicedatqc%mstp(jobs), & 1639 & ' ntyp = ', seaicedatqc%ntyp(jobs) 1640 ENDIF 1641 CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' ) 1642 1643 ENDIF 1644 1645 zlam = seaicedatqc%rlam(jobs) 1646 zphi = seaicedatqc%rphi(jobs) 1647 1648 ! Get weights to interpolate the model sea ice to the observation point 1649 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 1650 & zglam(:,:,iobs), zgphi(:,:,iobs), & 1651 & zmask(:,:,iobs), zweig, zobsmask ) 1652 1653 ! ... Interpolate the model sea ice to the observation point 1654 CALL obs_int_h2d( 1, 1, & 1655 & zweig, zseaicel(:,:,iobs), zext ) 1656 1657 seaicedatqc%rmod(jobs,1) = zext(1) 1658 1659 END DO 1660 1661 ! Deallocate the data for interpolation 1662 DEALLOCATE( & 1663 & igrdi, & 1664 & igrdj, & 1665 & zglam, & 1666 & zgphi, & 1667 & zmask, & 1668 & zseaicel & 1669 & ) 1670 1671 seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice 1672 1673 END SUBROUTINE obs_seaice_opt 1674 1675 SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 1676 & pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, & 1677 & ld_dailyav ) 1678 !!----------------------------------------------------------------------- 1679 !! 1680 !! *** ROUTINE obs_vel_opt *** 1681 !! 1682 !! ** Purpose : Compute the model counterpart of velocity profile 1683 !! data by interpolating from the model grid to the 1684 !! observation point. 1685 !! 1686 !! ** Method : Linearly interpolate zonal and meridional components of velocity 1687 !! to each observation point using the model values at the corners of 1688 !! the surrounding grid box. The model velocity components are on a 1689 !! staggered C- grid. 1690 !! 1691 !! For velocity data from the TAO array, the model equivalent is 1692 !! a daily mean velocity field. So, we first compute 1693 !! the mean, then interpolate only at the end of the day. 1694 !! 1695 !! ** Action : 1696 !! 1697 !! History : 1698 !! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles 1699 !! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles 1700 !!----------------------------------------------------------------------- 1701 1702 !! * Modules used 1703 USE obs_profiles_def ! Definition of storage space for profile obs. 1704 1705 IMPLICIT NONE 1706 1707 !! * Arguments 1708 TYPE(obs_prof), INTENT(INOUT) :: & 1709 & prodatqc ! Subset of profile data not failing screening 1710 INTEGER, INTENT(IN) :: kt ! Time step 1711 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1712 INTEGER, INTENT(IN) :: kpj 1713 INTEGER, INTENT(IN) :: kpk 1714 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1715 ! (kit000-1 = restart time) 1716 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 1717 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1718 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 1719 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 1720 & pun, & ! Model zonal component of velocity 1721 & pvn, & ! Model meridional component of velocity 1722 & pumask, & ! Land-sea mask 1723 & pvmask ! Land-sea mask 1724 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 1725 & pgdept ! Model array of depth levels 1726 LOGICAL, INTENT(IN) :: ld_dailyav 1727 1728 !! * Local declarations 1729 INTEGER :: ji 1730 INTEGER :: jj 1731 INTEGER :: jk 1732 INTEGER :: jobs 1733 INTEGER :: inrc 1734 INTEGER :: ipro 1735 INTEGER :: idayend 1736 INTEGER :: ista 1737 INTEGER :: iend 1738 INTEGER :: iobs 1739 INTEGER, DIMENSION(imaxavtypes) :: & 1740 & idailyavtypes 1741 REAL(KIND=wp) :: zlam 1742 REAL(KIND=wp) :: zphi 1743 REAL(KIND=wp) :: zdaystp 1744 REAL(KIND=wp), DIMENSION(kpk) :: & 1745 & zobsmasku, & 1746 & zobsmaskv, & 1747 & zobsmask, & 1748 & zobsk, & 1749 & zobs2k 1750 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 1751 & zweigu,zweigv 1752 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 1753 & zumask, zvmask, & 1754 & zintu, & 1755 & zintv, & 1756 & zinmu, & 1757 & zinmv 1758 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1759 & zglamu, zglamv, & 1760 & zgphiu, zgphiv 1761 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1762 & igrdiu, & 1763 & igrdju, & 1764 & igrdiv, & 1765 & igrdjv 1766 1767 !------------------------------------------------------------------------ 1768 ! Local initialization 1769 !------------------------------------------------------------------------ 1770 ! ... Record and data counters 1771 inrc = kt - kit000 + 2 1772 ipro = prodatqc%npstp(inrc) 1773 1774 ! Initialize daily mean for first timestep 1775 idayend = MOD( kt - kit000 + 1, kdaystp ) 1776 1777 ! Added kt == 0 test to catch restart case 1778 IF ( idayend == 1 .OR. kt == 0) THEN 1779 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 1780 prodatqc%vdmean(:,:,:,1) = 0.0 1781 prodatqc%vdmean(:,:,:,2) = 0.0 1782 ENDIF 1783 1784 ! Increment the zonal velocity field for computing daily mean 1785 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:) 1786 ! Increment the meridional velocity field for computing daily mean 1787 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:) 1788 1789 ! Compute the daily mean at the end of day 1790 zdaystp = 1.0 / REAL( kdaystp ) 1791 IF ( idayend == 0 ) THEN 1792 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp 1793 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp 1794 ENDIF 1795 1796 ! Get the data for interpolation 1797 ALLOCATE( & 1798 & igrdiu(2,2,ipro), & 1799 & igrdju(2,2,ipro), & 1800 & igrdiv(2,2,ipro), & 1801 & igrdjv(2,2,ipro), & 1802 & zglamu(2,2,ipro), zglamv(2,2,ipro), & 1803 & zgphiu(2,2,ipro), zgphiv(2,2,ipro), & 1804 & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), & 1805 & zintu(2,2,kpk,ipro), & 1806 & zintv(2,2,kpk,ipro) & 1807 & ) 1808 1809 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1810 iobs = jobs - prodatqc%nprofup 1811 igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1 1812 igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1 1813 igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1 1814 igrdju(1,2,iobs) = prodatqc%mj(jobs,1) 1815 igrdiu(2,1,iobs) = prodatqc%mi(jobs,1) 1816 igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1 1817 igrdiu(2,2,iobs) = prodatqc%mi(jobs,1) 1818 igrdju(2,2,iobs) = prodatqc%mj(jobs,1) 1819 igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1 1820 igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1 1821 igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1 1822 igrdjv(1,2,iobs) = prodatqc%mj(jobs,2) 1823 igrdiv(2,1,iobs) = prodatqc%mi(jobs,2) 1824 igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1 1825 igrdiv(2,2,iobs) = prodatqc%mi(jobs,2) 1826 igrdjv(2,2,iobs) = prodatqc%mj(jobs,2) 1827 END DO 1828 1829 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu ) 1830 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu ) 1831 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask ) 1832 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu ) 1833 1834 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv ) 1835 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv ) 1836 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask ) 1837 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv ) 1838 1839 ! At the end of the day also get interpolated means 1840 IF ( idayend == 0 ) THEN 1841 1842 ALLOCATE( & 1843 & zinmu(2,2,kpk,ipro), & 1844 & zinmv(2,2,kpk,ipro) & 1845 & ) 1846 1847 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, & 1848 & prodatqc%vdmean(:,:,:,1), zinmu ) 1849 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, & 1850 & prodatqc%vdmean(:,:,:,2), zinmv ) 1851 1852 ENDIF 1853 1854 ! loop over observations 1855 1856 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1857 1858 iobs = jobs - prodatqc%nprofup 1859 1860 IF ( kt /= prodatqc%mstp(jobs) ) THEN 1861 1862 IF(lwp) THEN 1863 WRITE(numout,*) 1864 WRITE(numout,*) ' E R R O R : Observation', & 1865 & ' time step is not consistent with the', & 1866 & ' model time step' 1867 WRITE(numout,*) ' =========' 1868 WRITE(numout,*) 1869 WRITE(numout,*) ' Record = ', jobs, & 1870 & ' kt = ', kt, & 1871 & ' mstp = ', prodatqc%mstp(jobs), & 1872 & ' ntyp = ', prodatqc%ntyp(jobs) 1873 ENDIF 1874 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 1875 ENDIF 1876 1877 zlam = prodatqc%rlam(jobs) 1878 zphi = prodatqc%rphi(jobs) 1879 1880 ! Initialize observation masks 1881 1882 zobsmasku(:) = 0.0 1883 zobsmaskv(:) = 0.0 1884 1885 ! Horizontal weights and vertical mask 1886 1887 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1888 1889 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1890 & zglamu(:,:,iobs), zgphiu(:,:,iobs), & 1891 & zumask(:,:,:,iobs), zweigu, zobsmasku ) 1892 1893 ENDIF 1894 1895 1896 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1897 1898 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1899 & zglamv(:,:,iobs), zgphiv(:,:,iobs), & 1900 & zvmask(:,:,:,iobs), zweigv, zobsmaskv ) 1901 1902 ENDIF 1903 1904 ! Ensure that the vertical mask on u and v are consistent. 1905 1906 zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) ) 1907 1908 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1909 1910 zobsk(:) = obfillflt 1911 1912 IF ( ld_dailyav ) THEN 1913 1914 IF ( idayend == 0 ) THEN 1915 1916 ! Daily averaged data 1917 1918 CALL obs_int_h2d( kpk, kpk, & 1919 & zweigu, zinmu(:,:,:,iobs), zobsk ) 1920 1921 1922 ELSE 1923 1924 CALL ctl_stop( ' A nonzero' // & 1925 & ' number of U profile data should' // & 1926 & ' only occur at the end of a given day' ) 1927 1928 ENDIF 1929 1930 ELSE 1931 1932 ! Point data 1933 1934 CALL obs_int_h2d( kpk, kpk, & 1935 & zweigu, zintu(:,:,:,iobs), zobsk ) 1936 1937 ENDIF 1938 1939 !------------------------------------------------------------- 1940 ! Compute vertical second-derivative of the interpolating 1941 ! polynomial at obs points 1942 !------------------------------------------------------------- 1943 1944 IF ( k1dint == 1 ) THEN 1945 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 1946 & pgdept, zobsmask ) 1947 ENDIF 1948 1949 !----------------------------------------------------------------- 1950 ! Vertical interpolation to the observation point 1951 !----------------------------------------------------------------- 1952 ista = prodatqc%npvsta(jobs,1) 1953 iend = prodatqc%npvend(jobs,1) 1954 CALL obs_int_z1d( kpk, & 1955 & prodatqc%var(1)%mvk(ista:iend), & 1956 & k1dint, iend - ista + 1, & 1957 & prodatqc%var(1)%vdep(ista:iend), & 1958 & zobsk, zobs2k, & 1959 & prodatqc%var(1)%vmod(ista:iend), & 1960 & pgdept, zobsmask ) 1961 1962 ENDIF 1963 1964 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1965 1966 zobsk(:) = obfillflt 1967 1968 IF ( ld_dailyav ) THEN 1969 1970 IF ( idayend == 0 ) THEN 1971 1972 ! Daily averaged data 1973 1974 CALL obs_int_h2d( kpk, kpk, & 1975 & zweigv, zinmv(:,:,:,iobs), zobsk ) 1976 1977 ELSE 1978 1979 CALL ctl_stop( ' A nonzero' // & 1980 & ' number of V profile data should' // & 1981 & ' only occur at the end of a given day' ) 1982 1983 ENDIF 1984 1985 ELSE 1986 1987 ! Point data 1988 1989 CALL obs_int_h2d( kpk, kpk, & 1990 & zweigv, zintv(:,:,:,iobs), zobsk ) 1991 1992 ENDIF 1993 1994 1995 !------------------------------------------------------------- 1996 ! Compute vertical second-derivative of the interpolating 1997 ! polynomial at obs points 1998 !------------------------------------------------------------- 1999 2000 IF ( k1dint == 1 ) THEN 2001 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 2002 & pgdept, zobsmask ) 2003 ENDIF 2004 2005 !---------------------------------------------------------------- 2006 ! Vertical interpolation to the observation point 2007 !---------------------------------------------------------------- 2008 ista = prodatqc%npvsta(jobs,2) 2009 iend = prodatqc%npvend(jobs,2) 2010 CALL obs_int_z1d( kpk, & 2011 & prodatqc%var(2)%mvk(ista:iend),& 2012 & k1dint, iend - ista + 1, & 2013 & prodatqc%var(2)%vdep(ista:iend),& 2014 & zobsk, zobs2k, & 2015 & prodatqc%var(2)%vmod(ista:iend),& 2016 & pgdept, zobsmask ) 2017 2018 ENDIF 2019 2020 END DO 2021 2022 ! Deallocate the data for interpolation 2023 DEALLOCATE( & 2024 & igrdiu, & 2025 & igrdju, & 2026 & igrdiv, & 2027 & igrdjv, & 2028 & zglamu, zglamv, & 2029 & zgphiu, zgphiv, & 2030 & zumask, zvmask, & 2031 & zintu, & 2032 & zintv & 2033 & ) 2034 ! At the end of the day also get interpolated means 2035 IF ( idayend == 0 ) THEN 2036 DEALLOCATE( & 2037 & zinmu, & 2038 & zinmv & 2039 & ) 2040 ENDIF 2041 2042 prodatqc%nprofup = prodatqc%nprofup + ipro 2043 2044 END SUBROUTINE obs_vel_opt 1341 1342 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1343 1344 END SUBROUTINE obs_surf_opt 2045 1345 2046 1346 END MODULE obs_oper 2047
Note: See TracChangeset
for help on using the changeset viewer.
