- Timestamp:
- 2015-12-16T16:44:35+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r4245 r6069 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 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and 12 !! salinity observations from profiles in generalised 13 !! vertical coordinates 22 14 !!---------------------------------------------------------------------- 23 15 24 !! * Modules used 16 !! * Modules used 25 17 USE par_kind, ONLY : & ! Precision variables 26 18 & wp 27 19 USE in_out_manager ! I/O manager 28 20 USE obs_inter_sup ! Interpolation support 29 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs ervationpt21 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 30 22 & obs_int_h2d, & 31 23 & obs_int_h2d_init 32 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs ervationpt24 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 33 25 & obs_int_z1d, & 34 26 & obs_int_z1d_spl … … 37 29 USE dom_oce, ONLY : & 38 30 & glamt, glamu, glamv, & 39 & gphit, gphiu, gphiv 31 & gphit, gphiu, gphiv, & 32 #if defined key_vvl 33 & gdept_n 34 #else 35 & gdept_0 36 #endif 40 37 USE lib_mpp, ONLY : & 41 38 & ctl_warn, ctl_stop 39 USE obs_grid, ONLY : & 40 & obs_level_search 41 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 42 & sbc_dcy, nday_qsr 42 43 43 44 IMPLICIT NONE … … 46 47 PRIVATE 47 48 48 PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations 49 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 50 & obs_sst_opt, & ! Compute the model counterpart of SST observations 51 & obs_sss_opt, & ! Compute the model counterpart of SSS observations 52 & obs_seaice_opt, & 53 & obs_vel_opt ! Compute the model counterpart of velocity profile data 54 55 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 49 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 50 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations 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 56 55 57 56 !!---------------------------------------------------------------------- … … 61 60 !!---------------------------------------------------------------------- 62 61 62 !! * Substitutions 63 # include "domzgr_substitute.h90" 63 64 CONTAINS 64 65 65 SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 66 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 67 & 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 68 72 !!----------------------------------------------------------------------- 69 73 !! … … 78 82 !! 79 83 !! First, a vertical profile of horizontally interpolated model 80 !! now temperatures is computed at the obs (lon, lat) point.84 !! now values is computed at the obs (lon, lat) point. 81 85 !! Several horizontal interpolation schemes are available: 82 86 !! - distance-weighted (great circle) (k2dint = 0) … … 86 90 !! - polynomial (quadrilateral grid) (k2dint = 4) 87 91 !! 88 !! Next, the vertical temperatureprofile is interpolated to the92 !! Next, the vertical profile is interpolated to the 89 93 !! data depth points. Two vertical interpolation schemes are 90 94 !! available: … … 96 100 !! routine. 97 101 !! 98 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is102 !! If the logical is switched on, the model equivalent is 99 103 !! a daily mean model temperature field. So, we first compute 100 104 !! the mean, then interpolate only at the end of the day. 101 105 !! 102 !! Note: thein situ temperature observations must be converted106 !! Note: in situ temperature observations must be converted 103 107 !! to potential temperature (the model variable) prior to 104 108 !! assimilation. 105 !!??????????????????????????????????????????????????????????????106 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???107 !!??????????????????????????????????????????????????????????????108 109 !! 109 110 !! ** Action : … … 115 116 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 116 117 !! ! 07-03 (K. Mogensen) General handling of profiles 118 !! ! 15-02 (M. Martin) Combined routine for all profile types 117 119 !!----------------------------------------------------------------------- 118 120 119 121 !! * Modules used 120 122 USE obs_profiles_def ! Definition of storage space for profile obs. … … 123 125 124 126 !! * Arguments 125 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 126 INTEGER, INTENT(IN) :: kt ! Time step 127 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 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 128 131 INTEGER, INTENT(IN) :: kpj 129 132 INTEGER, INTENT(IN) :: kpk 130 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step131 132 INTEGER, INTENT(IN) :: k1dint 133 INTEGER, INTENT(IN) :: k2dint 134 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per 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 135 138 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 136 & ptn, & ! Model temperature field 137 & psn, & ! Model salinity field 138 & ptmask ! Land-sea mask 139 & 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 139 148 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 140 & pgdept ! Model array of depth levels149 & pgdept ! Model array of depth levels 141 150 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 142 & kdailyavtypes! Types for daily averages 151 & kdailyavtypes ! Types for daily averages 152 143 153 !! * Local declarations 144 154 INTEGER :: ji … … 154 164 INTEGER, DIMENSION(imaxavtypes) :: & 155 165 & idailyavtypes 166 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 167 & igrdi1, & 168 & igrdi2, & 169 & igrdj1, & 170 & igrdj2 156 171 REAL(KIND=wp) :: zlam 157 172 REAL(KIND=wp) :: zphi 158 173 REAL(KIND=wp) :: zdaystp 159 174 REAL(KIND=wp), DIMENSION(kpk) :: & 160 & zobsmask, & 175 & zobsmask1, & 176 & zobsmask2, & 161 177 & zobsk, & 162 178 & zobs2k 163 179 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 164 & zweig 180 & zweig1, & 181 & zweig2 165 182 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 166 & zmask, & 167 & zintt, & 168 & zints, & 169 & zinmt, & 170 & zinms 183 & zmask1, & 184 & zmask2, & 185 & zint1, & 186 & zint2, & 187 & zinm1, & 188 & zinm2 171 189 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 172 & zglam , &173 & zg phi174 INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::&175 & igrdi, &176 & igrdj190 & zglam1, & 191 & zglam2, & 192 & zgphi1, & 193 & zgphi2 194 LOGICAL :: ld_dailyav 177 195 178 196 !------------------------------------------------------------------------ 179 197 ! Local initialization 180 198 !------------------------------------------------------------------------ 181 ! ...Record and data counters199 ! Record and data counters 182 200 inrc = kt - kit000 + 2 183 201 ipro = prodatqc%npstp(inrc) 184 202 185 203 ! Daily average types 204 ld_dailyav = .FALSE. 186 205 IF ( PRESENT(kdailyavtypes) ) THEN 187 206 idailyavtypes(:) = kdailyavtypes(:) 207 IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 188 208 ELSE 189 209 idailyavtypes(:) = -1 190 210 ENDIF 191 211 192 ! Initialize daily mean for first timestep 212 ! Daily means are calculated for values over timesteps: 213 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 193 214 idayend = MOD( kt - kit000 + 1, kdaystp ) 194 215 195 ! Added kt == 0 test to catch restart case 196 IF ( idayend == 1 .OR. kt == 0) THEN 197 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 216 IF ( ld_dailyav ) THEN 217 218 ! Initialize daily mean for first timestep of the day 219 IF ( idayend == 1 .OR. kt == 0 ) THEN 220 DO jk = 1, jpk 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 prodatqc%vdmean(ji,jj,jk,1) = 0.0 224 prodatqc%vdmean(ji,jj,jk,2) = 0.0 225 END DO 226 END DO 227 END DO 228 ENDIF 229 198 230 DO jk = 1, jpk 199 231 DO jj = 1, jpj 200 232 DO ji = 1, jpi 201 prodatqc%vdmean(ji,jj,jk,1) = 0.0 202 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) 203 239 END DO 204 240 END DO 205 241 END DO 206 ENDIF 207 208 DO jk = 1, jpk 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 ! Increment the temperature field for computing daily mean 212 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 213 & + ptn(ji,jj,jk) 214 ! Increment the salinity field for computing daily mean 215 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 216 & + psn(ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 ! Compute the daily mean at the end of day 222 zdaystp = 1.0 / REAL( kdaystp ) 223 IF ( idayend == 0 ) THEN 224 DO jk = 1, jpk 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 228 & * zdaystp 229 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 230 & * zdaystp 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 231 256 END DO 232 257 END DO 233 END DO 258 ENDIF 259 234 260 ENDIF 235 261 236 262 ! Get the data for interpolation 237 263 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) & 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) & 245 276 & ) 246 277 247 278 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 248 279 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) 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) 257 296 END DO 258 297 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 ) 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 ) 264 307 265 308 ! At the end of the day also get interpolated means 266 IF ( idayend == 0 ) THEN309 IF ( ld_dailyav .AND. idayend == 0 ) THEN 267 310 268 311 ALLOCATE( & 269 & zinm t(2,2,kpk,ipro), &270 & zinm s(2,2,kpk,ipro) &312 & zinm1(2,2,kpk,ipro), & 313 & zinm2(2,2,kpk,ipro) & 271 314 & ) 272 315 273 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &274 & prodatqc%vdmean(:,:,:,1), zinm t)275 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &276 & prodatqc%vdmean(:,:,:,2), zinm s)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 ) 277 320 278 321 ENDIF … … 283 326 284 327 IF ( kt /= prodatqc%mstp(jobs) ) THEN 285 328 286 329 IF(lwp) THEN 287 330 WRITE(numout,*) … … 298 341 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 299 342 ENDIF 300 343 301 344 zlam = prodatqc%rlam(jobs) 302 345 zphi = prodatqc%rphi(jobs) 303 346 304 347 ! Horizontal weights and vertical mask 305 348 306 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 307 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 349 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 308 350 309 351 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 310 & zglam (:,:,iobs), zgphi(:,:,iobs), &311 & zmask (:,:,:,iobs), zweig, zobsmask)352 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 353 & zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 312 354 313 355 ENDIF 314 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 315 365 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 316 366 317 367 zobsk(:) = obfillflt 318 368 319 369 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 320 370 321 371 IF ( idayend == 0 ) THEN 322 323 ! Daily averaged moored buoy (MRB) data 324 372 ! Daily averaged data 325 373 CALL obs_int_h2d( kpk, kpk, & 326 & zweig, zinmt(:,:,:,iobs), zobsk ) 327 328 329 ELSE 330 331 CALL ctl_stop( ' A nonzero' // & 332 & ' number of profile T BUOY data should' // & 333 & ' only occur at the end of a given day' ) 374 & zweig1, zinm1(:,:,:,iobs), zobsk ) 334 375 335 376 ENDIF 336 377 337 378 ELSE 338 379 339 380 ! Point data 340 341 381 CALL obs_int_h2d( kpk, kpk, & 342 & zweig , zintt(:,:,:,iobs), zobsk )382 & zweig1, zint1(:,:,:,iobs), zobsk ) 343 383 344 384 ENDIF … … 348 388 ! polynomial at obs points 349 389 !------------------------------------------------------------- 350 390 351 391 IF ( k1dint == 1 ) THEN 352 392 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 353 & pgdept, zobsmask )393 & pgdept, zobsmask1 ) 354 394 ENDIF 355 395 356 396 !----------------------------------------------------------------- 357 397 ! Vertical interpolation to the observation point … … 365 405 & zobsk, zobs2k, & 366 406 & prodatqc%var(1)%vmod(ista:iend), & 367 & pgdept, zobsmask )407 & pgdept, zobsmask1 ) 368 408 369 409 ENDIF … … 377 417 IF ( idayend == 0 ) THEN 378 418 379 ! Daily averaged moored buoy (MRB) data 380 419 ! Daily averaged data 381 420 CALL obs_int_h2d( kpk, kpk, & 382 & zweig, zinms(:,:,:,iobs), zobsk ) 383 384 ELSE 385 386 CALL ctl_stop( ' A nonzero' // & 387 & ' number of profile S BUOY data should' // & 388 & ' only occur at the end of a given day' ) 421 & zweig2, zinm2(:,:,:,iobs), zobsk ) 389 422 390 423 ENDIF 391 424 392 425 ELSE 393 426 394 427 ! Point data 395 396 428 CALL obs_int_h2d( kpk, kpk, & 397 & zweig , zints(:,:,:,iobs), zobsk )429 & zweig2, zint2(:,:,:,iobs), zobsk ) 398 430 399 431 ENDIF … … 404 436 ! polynomial at obs points 405 437 !------------------------------------------------------------- 406 438 407 439 IF ( k1dint == 1 ) THEN 408 440 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 409 & pgdept, zobsmask )441 & pgdept, zobsmask2 ) 410 442 ENDIF 411 443 412 444 !---------------------------------------------------------------- 413 445 ! Vertical interpolation to the observation point … … 421 453 & zobsk, zobs2k, & 422 454 & prodatqc%var(2)%vmod(ista:iend),& 423 & pgdept, zobsmask )455 & pgdept, zobsmask2 ) 424 456 425 457 ENDIF 426 458 427 459 END DO 428 460 429 461 ! Deallocate the data for interpolation 430 462 DEALLOCATE( & 431 & igrdi, & 432 & igrdj, & 433 & zglam, & 434 & zgphi, & 435 & zmask, & 436 & zintt, & 437 & 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 & 438 475 & ) 476 439 477 ! At the end of the day also get interpolated means 440 IF ( idayend == 0 ) THEN478 IF ( ld_dailyav .AND. idayend == 0 ) THEN 441 479 DEALLOCATE( & 442 & zinm t, &443 & zinm s&480 & zinm1, & 481 & zinm2 & 444 482 & ) 445 483 ENDIF 446 484 447 485 prodatqc%nprofup = prodatqc%nprofup + ipro 486 487 END SUBROUTINE obs_prof_opt 488 489 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 490 & ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 491 & kdailyavtypes ) 492 !!----------------------------------------------------------------------- 493 !! 494 !! *** ROUTINE obs_pro_opt *** 495 !! 496 !! ** Purpose : Compute the model counterpart of profiles 497 !! data by interpolating from the model grid to the 498 !! observation point. Generalised vertical coordinate version 499 !! 500 !! ** Method : Linearly interpolate to each observation point using 501 !! the model values at the corners of the surrounding grid box. 502 !! 503 !! First, model values on the model grid are interpolated vertically to the 504 !! Depths of the profile observations. Two vertical interpolation schemes are 505 !! available: 506 !! - linear (k1dint = 0) 507 !! - Cubic spline (k1dint = 1) 508 !! 509 !! 510 !! Secondly the interpolated values are interpolated horizontally to the 511 !! obs (lon, lat) point. 512 !! Several horizontal interpolation schemes are available: 513 !! - distance-weighted (great circle) (k2dint = 0) 514 !! - distance-weighted (small angle) (k2dint = 1) 515 !! - bilinear (geographical grid) (k2dint = 2) 516 !! - bilinear (quadrilateral grid) (k2dint = 3) 517 !! - polynomial (quadrilateral grid) (k2dint = 4) 518 !! 519 !! For the cubic spline the 2nd derivative of the interpolating 520 !! polynomial is computed before entering the vertical interpolation 521 !! routine. 522 !! 523 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is 524 !! a daily mean model temperature field. So, we first compute 525 !! the mean, then interpolate only at the end of the day. 526 !! 527 !! This is the procedure to be used with generalised vertical model 528 !! coordinates (ie s-coordinates. It is ~4x slower than the equivalent 529 !! horizontal then vertical interpolation algorithm, but can deal with situations 530 !! where the model levels are not flat. 531 !! ONLY PERFORMED if ln_sco=.TRUE. 532 !! 533 !! Note: the in situ temperature observations must be converted 534 !! to potential temperature (the model variable) prior to 535 !! assimilation. 536 !!?????????????????????????????????????????????????????????????? 537 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 538 !!?????????????????????????????????????????????????????????????? 539 !! 540 !! ** Action : 541 !! 542 !! History : 543 !! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 544 !! vertical coordinates 545 !!----------------------------------------------------------------------- 546 547 !! * Modules used 548 USE obs_profiles_def ! Definition of storage space for profile obs. 549 USE dom_oce, ONLY : & 550 #if defined key_vvl 551 & gdepw_n 552 #else 553 & gdepw_0 554 #endif 555 556 IMPLICIT NONE 557 558 !! * Arguments 559 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 560 INTEGER, INTENT(IN) :: kt ! Time step 561 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 562 INTEGER, INTENT(IN) :: kpj 563 INTEGER, INTENT(IN) :: kpk 564 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 565 ! (kit000-1 = restart time) 566 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 567 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 568 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 569 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 570 & ptn, & ! Model temperature field 571 & psn, & ! Model salinity field 572 & ptmask ! Land-sea mask 573 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 574 & pgdept, & ! Model array of depth T levels 575 & pgdepw ! Model array of depth W levels 576 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 577 & kdailyavtypes ! Types for daily averages 448 578 449 END SUBROUTINE obs_pro_opt 450 451 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 & psshn, psshmask, k2dint ) 579 !! * Local declarations 580 INTEGER :: ji 581 INTEGER :: jj 582 INTEGER :: jk 583 INTEGER :: iico, ijco 584 INTEGER :: jobs 585 INTEGER :: inrc 586 INTEGER :: ipro 587 INTEGER :: idayend 588 INTEGER :: ista 589 INTEGER :: iend 590 INTEGER :: iobs 591 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 592 INTEGER, DIMENSION(imaxavtypes) :: & 593 & idailyavtypes 594 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 595 & igrdi, & 596 & igrdj 597 INTEGER :: & 598 & inum_obs 599 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 600 REAL(KIND=wp) :: zlam 601 REAL(KIND=wp) :: zphi 602 REAL(KIND=wp) :: zdaystp 603 REAL(KIND=wp), DIMENSION(kpk) :: & 604 & zobsmask, & 605 & zobsk, & 606 & zobs2k 607 REAL(KIND=wp), DIMENSION(2,2,1) :: & 608 & zweig, & 609 & l_zweig 610 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 611 & zmask, & 612 & zintt, & 613 & zints, & 614 & zinmt, & 615 & zgdept,& 616 & zgdepw,& 617 & zinms 618 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 619 & zglam, & 620 & zgphi 621 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 622 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 623 624 !------------------------------------------------------------------------ 625 ! Local initialization 626 !------------------------------------------------------------------------ 627 ! ... Record and data counters 628 inrc = kt - kit000 + 2 629 ipro = prodatqc%npstp(inrc) 630 631 ! Daily average types 632 IF ( PRESENT(kdailyavtypes) ) THEN 633 idailyavtypes(:) = kdailyavtypes(:) 634 ELSE 635 idailyavtypes(:) = -1 636 ENDIF 637 638 ! Initialize daily mean for first time-step 639 idayend = MOD( kt - kit000 + 1, kdaystp ) 640 641 ! Added kt == 0 test to catch restart case 642 IF ( idayend == 1 .OR. kt == 0) THEN 643 644 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 645 DO jk = 1, jpk 646 DO jj = 1, jpj 647 DO ji = 1, jpi 648 prodatqc%vdmean(ji,jj,jk,1) = 0.0 649 prodatqc%vdmean(ji,jj,jk,2) = 0.0 650 END DO 651 END DO 652 END DO 653 654 ENDIF 655 656 DO jk = 1, jpk 657 DO jj = 1, jpj 658 DO ji = 1, jpi 659 ! Increment the temperature field for computing daily mean 660 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 661 & + ptn(ji,jj,jk) 662 ! Increment the salinity field for computing daily mean 663 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 664 & + psn(ji,jj,jk) 665 END DO 666 END DO 667 END DO 668 669 ! Compute the daily mean at the end of day 670 zdaystp = 1.0 / REAL( kdaystp ) 671 IF ( idayend == 0 ) THEN 672 DO jk = 1, jpk 673 DO jj = 1, jpj 674 DO ji = 1, jpi 675 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 676 & * zdaystp 677 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 678 & * zdaystp 679 END DO 680 END DO 681 END DO 682 ENDIF 683 684 ! Get the data for interpolation 685 ALLOCATE( & 686 & igrdi(2,2,ipro), & 687 & igrdj(2,2,ipro), & 688 & zglam(2,2,ipro), & 689 & zgphi(2,2,ipro), & 690 & zmask(2,2,kpk,ipro), & 691 & zintt(2,2,kpk,ipro), & 692 & zints(2,2,kpk,ipro), & 693 & zgdept(2,2,kpk,ipro), & 694 & zgdepw(2,2,kpk,ipro) & 695 & ) 696 697 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 698 iobs = jobs - prodatqc%nprofup 699 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 700 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 701 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 702 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 703 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 704 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 705 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 706 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 707 END DO 708 709 ! Initialise depth arrays 710 zgdept = 0.0 711 zgdepw = 0.0 712 713 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam ) 714 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi ) 715 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask ) 716 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn, zintt ) 717 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn, zints ) 718 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), & 719 & zgdept ) 720 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), & 721 & zgdepw ) 722 723 ! At the end of the day also get interpolated means 724 IF ( idayend == 0 ) THEN 725 726 ALLOCATE( & 727 & zinmt(2,2,kpk,ipro), & 728 & zinms(2,2,kpk,ipro) & 729 & ) 730 731 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 732 & prodatqc%vdmean(:,:,:,1), zinmt ) 733 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 734 & prodatqc%vdmean(:,:,:,2), zinms ) 735 736 ENDIF 737 738 ! Return if no observations to process 739 ! Has to be done after comm commands to ensure processors 740 ! stay in sync 741 IF ( ipro == 0 ) RETURN 742 743 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 744 745 iobs = jobs - prodatqc%nprofup 746 747 IF ( kt /= prodatqc%mstp(jobs) ) THEN 748 749 IF(lwp) THEN 750 WRITE(numout,*) 751 WRITE(numout,*) ' E R R O R : Observation', & 752 & ' time step is not consistent with the', & 753 & ' model time step' 754 WRITE(numout,*) ' =========' 755 WRITE(numout,*) 756 WRITE(numout,*) ' Record = ', jobs, & 757 & ' kt = ', kt, & 758 & ' mstp = ', prodatqc%mstp(jobs), & 759 & ' ntyp = ', prodatqc%ntyp(jobs) 760 ENDIF 761 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 762 ENDIF 763 764 zlam = prodatqc%rlam(jobs) 765 zphi = prodatqc%rphi(jobs) 766 767 ! Horizontal weights 768 ! Only calculated once, for both T and S. 769 ! Masked values are calculated later. 770 771 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 772 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 773 774 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 775 & zglam(:,:,iobs), zgphi(:,:,iobs), & 776 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 777 778 ENDIF 779 780 ! IF zmsk_1 = 0; then ob is on land 781 IF (zmsk_1(1) < 0.1) THEN 782 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 783 784 ELSE 785 786 ! Temperature 787 788 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 789 790 zobsk(:) = obfillflt 791 792 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 793 794 IF ( idayend == 0 ) THEN 795 796 ! Daily averaged moored buoy (MRB) data 797 798 ! vertically interpolate all 4 corners 799 ista = prodatqc%npvsta(jobs,1) 800 iend = prodatqc%npvend(jobs,1) 801 inum_obs = iend - ista + 1 802 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 803 804 DO iin=1,2 805 DO ijn=1,2 806 807 808 809 IF ( k1dint == 1 ) THEN 810 CALL obs_int_z1d_spl( kpk, & 811 & zinmt(iin,ijn,:,iobs), & 812 & zobs2k, zgdept(iin,ijn,:,iobs), & 813 & zmask(iin,ijn,:,iobs)) 814 ENDIF 815 816 CALL obs_level_search(kpk, & 817 & zgdept(iin,ijn,:,iobs), & 818 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 819 & iv_indic) 820 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 821 & prodatqc%var(1)%vdep(ista:iend), & 822 & zinmt(iin,ijn,:,iobs), & 823 & zobs2k, interp_corner(iin,ijn,:), & 824 & zgdept(iin,ijn,:,iobs), & 825 & zmask(iin,ijn,:,iobs)) 826 827 ENDDO 828 ENDDO 829 830 831 ELSE 832 833 CALL ctl_stop( ' A nonzero' // & 834 & ' number of profile T BUOY data should' // & 835 & ' only occur at the end of a given day' ) 836 837 ENDIF 838 839 ELSE 840 841 ! Point data 842 843 ! vertically interpolate all 4 corners 844 ista = prodatqc%npvsta(jobs,1) 845 iend = prodatqc%npvend(jobs,1) 846 inum_obs = iend - ista + 1 847 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 848 DO iin=1,2 849 DO ijn=1,2 850 851 852 IF ( k1dint == 1 ) THEN 853 CALL obs_int_z1d_spl( kpk, & 854 & zintt(iin,ijn,:,iobs),& 855 & zobs2k, zgdept(iin,ijn,:,iobs), & 856 & zmask(iin,ijn,:,iobs)) 857 858 ENDIF 859 860 CALL obs_level_search(kpk, & 861 & zgdept(iin,ijn,:,iobs),& 862 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 863 & iv_indic) 864 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 865 & prodatqc%var(1)%vdep(ista:iend), & 866 & zintt(iin,ijn,:,iobs), & 867 & zobs2k,interp_corner(iin,ijn,:), & 868 & zgdept(iin,ijn,:,iobs), & 869 & zmask(iin,ijn,:,iobs) ) 870 871 ENDDO 872 ENDDO 873 874 ENDIF 875 876 !------------------------------------------------------------- 877 ! Compute the horizontal interpolation for every profile level 878 !------------------------------------------------------------- 879 880 DO ikn=1,inum_obs 881 iend=ista+ikn-1 882 883 l_zweig(:,:,1) = 0._wp 884 885 ! This code forces the horizontal weights to be 886 ! zero IF the observation is below the bottom of the 887 ! corners of the interpolation nodes, Or if it is in 888 ! the mask. This is important for observations are near 889 ! steep bathymetry 890 DO iin=1,2 891 DO ijn=1,2 892 893 depth_loop1: DO ik=kpk,2,-1 894 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 895 896 l_zweig(iin,ijn,1) = & 897 & zweig(iin,ijn,1) * & 898 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 899 & - prodatqc%var(1)%vdep(iend)),0._wp) 900 901 EXIT depth_loop1 902 ENDIF 903 ENDDO depth_loop1 904 905 ENDDO 906 ENDDO 907 908 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 909 & prodatqc%var(1)%vmod(iend:iend) ) 910 911 ENDDO 912 913 914 DEALLOCATE(interp_corner,iv_indic) 915 916 ENDIF 917 918 919 ! Salinity 920 921 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 922 923 zobsk(:) = obfillflt 924 925 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 926 927 IF ( idayend == 0 ) THEN 928 929 ! Daily averaged moored buoy (MRB) data 930 931 ! vertically interpolate all 4 corners 932 ista = prodatqc%npvsta(iobs,2) 933 iend = prodatqc%npvend(iobs,2) 934 inum_obs = iend - ista + 1 935 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 936 937 DO iin=1,2 938 DO ijn=1,2 939 940 941 942 IF ( k1dint == 1 ) THEN 943 CALL obs_int_z1d_spl( kpk, & 944 & zinms(iin,ijn,:,iobs), & 945 & zobs2k, zgdept(iin,ijn,:,iobs), & 946 & zmask(iin,ijn,:,iobs)) 947 ENDIF 948 949 CALL obs_level_search(kpk, & 950 & zgdept(iin,ijn,:,iobs), & 951 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 952 & iv_indic) 953 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 954 & prodatqc%var(2)%vdep(ista:iend), & 955 & zinms(iin,ijn,:,iobs), & 956 & zobs2k, interp_corner(iin,ijn,:), & 957 & zgdept(iin,ijn,:,iobs), & 958 & zmask(iin,ijn,:,iobs)) 959 960 ENDDO 961 ENDDO 962 963 964 ELSE 965 966 CALL ctl_stop( ' A nonzero' // & 967 & ' number of profile T BUOY data should' // & 968 & ' only occur at the end of a given day' ) 969 970 ENDIF 971 972 ELSE 973 974 ! Point data 975 976 ! vertically interpolate all 4 corners 977 ista = prodatqc%npvsta(jobs,2) 978 iend = prodatqc%npvend(jobs,2) 979 inum_obs = iend - ista + 1 980 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 981 982 DO iin=1,2 983 DO ijn=1,2 984 985 986 IF ( k1dint == 1 ) THEN 987 CALL obs_int_z1d_spl( kpk, & 988 & zints(iin,ijn,:,iobs),& 989 & zobs2k, zgdept(iin,ijn,:,iobs), & 990 & zmask(iin,ijn,:,iobs)) 991 992 ENDIF 993 994 CALL obs_level_search(kpk, & 995 & zgdept(iin,ijn,:,iobs),& 996 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 997 & iv_indic) 998 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 999 & prodatqc%var(2)%vdep(ista:iend), & 1000 & zints(iin,ijn,:,iobs), & 1001 & zobs2k,interp_corner(iin,ijn,:), & 1002 & zgdept(iin,ijn,:,iobs), & 1003 & zmask(iin,ijn,:,iobs) ) 1004 1005 ENDDO 1006 ENDDO 1007 1008 ENDIF 1009 1010 !------------------------------------------------------------- 1011 ! Compute the horizontal interpolation for every profile level 1012 !------------------------------------------------------------- 1013 1014 DO ikn=1,inum_obs 1015 iend=ista+ikn-1 1016 1017 l_zweig(:,:,1) = 0._wp 1018 1019 ! This code forces the horizontal weights to be 1020 ! zero IF the observation is below the bottom of the 1021 ! corners of the interpolation nodes, Or if it is in 1022 ! the mask. This is important for observations are near 1023 ! steep bathymetry 1024 DO iin=1,2 1025 DO ijn=1,2 1026 1027 depth_loop2: DO ik=kpk,2,-1 1028 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 1029 1030 l_zweig(iin,ijn,1) = & 1031 & zweig(iin,ijn,1) * & 1032 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 1033 & - prodatqc%var(2)%vdep(iend)),0._wp) 1034 1035 EXIT depth_loop2 1036 ENDIF 1037 ENDDO depth_loop2 1038 1039 ENDDO 1040 ENDDO 1041 1042 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1043 & prodatqc%var(2)%vmod(iend:iend) ) 1044 1045 ENDDO 1046 1047 1048 DEALLOCATE(interp_corner,iv_indic) 1049 1050 ENDIF 1051 1052 ENDIF 1053 1054 END DO 1055 1056 ! Deallocate the data for interpolation 1057 DEALLOCATE( & 1058 & igrdi, & 1059 & igrdj, & 1060 & zglam, & 1061 & zgphi, & 1062 & zmask, & 1063 & zintt, & 1064 & zints, & 1065 & zgdept,& 1066 & zgdepw & 1067 & ) 1068 ! At the end of the day also get interpolated means 1069 IF ( idayend == 0 ) THEN 1070 DEALLOCATE( & 1071 & zinmt, & 1072 & zinms & 1073 & ) 1074 ENDIF 1075 1076 prodatqc%nprofup = prodatqc%nprofup + ipro 1077 1078 END SUBROUTINE obs_pro_sco_opt 1079 1080 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 1081 & kit000, kdaystp, psurf, psurfmask, & 1082 & k2dint, ldnightav ) 1083 453 1084 !!----------------------------------------------------------------------- 454 1085 !! 455 !! *** ROUTINE obs_s la_opt ***456 !! 457 !! ** Purpose : Compute the model counterpart of s ea level anomaly1086 !! *** ROUTINE obs_surf_opt *** 1087 !! 1088 !! ** Purpose : Compute the model counterpart of surface 458 1089 !! data by interpolating from the model grid to the 459 1090 !! observation point. … … 462 1093 !! the model values at the corners of the surrounding grid box. 463 1094 !! 464 !! 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. 465 1096 !! 466 1097 !! Several horizontal interpolation schemes are available: … … 470 1101 !! - bilinear (quadrilateral grid) (k2dint = 3) 471 1102 !! - polynomial (quadrilateral grid) (k2dint = 4) 472 !! 473 !! The sea level anomaly at the observation points is then computed 474 !! by removing a mean dynamic topography (defined at the obs. point). 1103 !! 475 1104 !! 476 1105 !! ** Action : … … 478 1107 !! History : 479 1108 !! ! 07-03 (A. Weaver) 1109 !! ! 15-02 (M. Martin) Combined routine for surface types 480 1110 !!----------------------------------------------------------------------- 481 1111 482 1112 !! * Modules used 483 1113 USE obs_surf_def ! Definition of storage space for surface observations … … 486 1116 487 1117 !! * Arguments 488 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening 489 INTEGER, INTENT(IN) :: kt ! Time step 490 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 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 491 1122 INTEGER, INTENT(IN) :: kpj 492 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 493 ! (kit000-1 = restart time) 494 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 495 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 496 & psshn, & ! Model SSH field 497 & psshmask ! Land-sea mask 498 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 499 1132 !! * Local declarations 500 1133 INTEGER :: ji … … 502 1135 INTEGER :: jobs 503 1136 INTEGER :: inrc 504 INTEGER :: is la1137 INTEGER :: isurf 505 1138 INTEGER :: iobs 506 REAL(KIND=wp) :: zlam 507 REAL(KIND=wp) :: zphi 508 REAL(KIND=wp) :: zext(1), zobsmask(1) 509 REAL(kind=wp), DIMENSION(2,2,1) :: & 510 & zweig 511 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 512 & zmask, & 513 & zsshl, & 514 & zglam, & 515 & zgphi 1139 INTEGER :: idayend 516 1140 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 517 1141 & igrdi, & 518 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) 519 1162 520 1163 !------------------------------------------------------------------------ 521 1164 ! Local initialization 522 1165 !------------------------------------------------------------------------ 523 ! ...Record and data counters1166 ! Record and data counters 524 1167 inrc = kt - kit000 + 2 525 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 526 1231 527 1232 ! Get the data for interpolation 528 1233 529 1234 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) &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) & 536 1241 & ) 537 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)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) 548 1253 END DO 549 1254 550 CALL obs_int_comm_2d( 2, 2, is la, &1255 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 551 1256 & igrdi, igrdj, glamt, zglam ) 552 CALL obs_int_comm_2d( 2, 2, is la, &1257 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 553 1258 & 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 ) 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 558 1275 559 1276 ! Loop over observations 560 561 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 562 563 iobs = jobs - sladatqc%nsurfup 564 565 IF ( kt /= sladatqc%mstp(jobs) ) THEN 566 1277 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1278 1279 iobs = jobs - surfdataqc%nsurfup 1280 1281 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 1282 567 1283 IF(lwp) THEN 568 1284 WRITE(numout,*) … … 574 1290 WRITE(numout,*) ' Record = ', jobs, & 575 1291 & ' kt = ', kt, & 576 & ' mstp = ', s ladatqc%mstp(jobs), &577 & ' ntyp = ', s ladatqc%ntyp(jobs)1292 & ' mstp = ', surfdataqc%mstp(jobs), & 1293 & ' ntyp = ', surfdataqc%ntyp(jobs) 578 1294 ENDIF 579 CALL ctl_stop( 'obs_s la_opt', 'Inconsistent time' )580 1295 CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 1296 581 1297 ENDIF 582 583 zlam = s ladatqc%rlam(jobs)584 zphi = s ladatqc%rphi(jobs)585 586 ! 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 587 1303 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 588 1304 & zglam(:,:,iobs), zgphi(:,:,iobs), & 589 1305 & zmask(:,:,iobs), zweig, zobsmask ) 590 591 592 ! Interpolate the model SSH to the observation point 593 CALL obs_int_h2d( 1, 1, & 594 & zweig, zsshl(:,:,iobs), zext ) 595 596 sladatqc%rext(jobs,1) = zext(1) 597 ! ... Remove the MDT at the observation point 598 sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 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 599 1322 600 1323 END DO … … 607 1330 & zgphi, & 608 1331 & zmask, & 609 & zs shl&1332 & zsurf & 610 1333 & ) 611 1334 612 sladatqc%nsurfup = sladatqc%nsurfup + isla 613 614 END SUBROUTINE obs_sla_opt 615 616 SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 617 & psstn, psstmask, k2dint, ld_nightav ) 618 !!----------------------------------------------------------------------- 619 !! 620 !! *** ROUTINE obs_sst_opt *** 621 !! 622 !! ** Purpose : Compute the model counterpart of surface temperature 623 !! data by interpolating from the model grid to the 624 !! observation point. 625 !! 626 !! ** Method : Linearly interpolate to each observation point using 627 !! the model values at the corners of the surrounding grid box. 628 !! 629 !! The now model SST is first computed at the obs (lon, lat) point. 630 !! 631 !! Several horizontal interpolation schemes are available: 632 !! - distance-weighted (great circle) (k2dint = 0) 633 !! - distance-weighted (small angle) (k2dint = 1) 634 !! - bilinear (geographical grid) (k2dint = 2) 635 !! - bilinear (quadrilateral grid) (k2dint = 3) 636 !! - polynomial (quadrilateral grid) (k2dint = 4) 637 !! 638 !! 639 !! ** Action : 640 !! 641 !! History : 642 !! ! 07-07 (S. Ricci ) : Original 643 !! 644 !!----------------------------------------------------------------------- 645 646 !! * Modules used 647 USE obs_surf_def ! Definition of storage space for surface observations 648 USE sbcdcy 649 650 IMPLICIT NONE 651 652 !! * Arguments 653 TYPE(obs_surf), INTENT(INOUT) :: & 654 & sstdatqc ! Subset of surface data not failing screening 655 INTEGER, INTENT(IN) :: kt ! Time step 656 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 657 INTEGER, INTENT(IN) :: kpj 658 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 659 ! (kit000-1 = restart time) 660 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 661 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 662 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 663 & psstn, & ! Model SST field 664 & psstmask ! Land-sea mask 665 666 !! * Local declarations 667 INTEGER :: ji 668 INTEGER :: jj 669 INTEGER :: jobs 670 INTEGER :: inrc 671 INTEGER :: isst 672 INTEGER :: iobs 673 INTEGER :: idayend 674 REAL(KIND=wp) :: zlam 675 REAL(KIND=wp) :: zphi 676 REAL(KIND=wp) :: zext(1), zobsmask(1) 677 REAL(KIND=wp) :: zdaystp 678 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 679 & icount_sstnight, & 680 & imask_night 681 REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 682 & zintmp, & 683 & zouttmp, & 684 & zmeanday ! to compute model sst in region of 24h daylight (pole) 685 REAL(kind=wp), DIMENSION(2,2,1) :: & 686 & zweig 687 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 688 & zmask, & 689 & zsstl, & 690 & zsstm, & 691 & zglam, & 692 & zgphi 693 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 694 & igrdi, & 695 & igrdj 696 LOGICAL, INTENT(IN) :: ld_nightav 697 698 !----------------------------------------------------------------------- 699 ! Local initialization 700 !----------------------------------------------------------------------- 701 ! ... Record and data counters 702 inrc = kt - kit000 + 2 703 isst = sstdatqc%nsstp(inrc) 704 705 IF ( ld_nightav ) THEN 706 707 ! Initialize array for night mean 708 709 IF ( kt .EQ. 0 ) THEN 710 ALLOCATE ( icount_sstnight(kpi,kpj) ) 711 ALLOCATE ( imask_night(kpi,kpj) ) 712 ALLOCATE ( zintmp(kpi,kpj) ) 713 ALLOCATE ( zouttmp(kpi,kpj) ) 714 ALLOCATE ( zmeanday(kpi,kpj) ) 715 nday_qsr = -1 ! initialisation flag for nbc_dcy 716 ENDIF 717 718 ! Initialize daily mean for first timestep 719 idayend = MOD( kt - kit000 + 1, kdaystp ) 720 721 ! Added kt == 0 test to catch restart case 722 IF ( idayend == 1 .OR. kt == 0) THEN 723 IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 724 DO jj = 1, jpj 725 DO ji = 1, jpi 726 sstdatqc%vdmean(ji,jj) = 0.0 727 zmeanday(ji,jj) = 0.0 728 icount_sstnight(ji,jj) = 0 729 END DO 730 END DO 731 ENDIF 732 733 zintmp(:,:) = 0.0 734 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 735 imask_night(:,:) = INT( zouttmp(:,:) ) 736 737 DO jj = 1, jpj 738 DO ji = 1, jpi 739 ! Increment the temperature field for computing night mean and counter 740 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 741 & + psstn(ji,jj)*imask_night(ji,jj) 742 zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj) 743 icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 744 END DO 745 END DO 746 747 ! Compute the daily mean at the end of day 748 749 zdaystp = 1.0 / REAL( kdaystp ) 750 751 IF ( idayend == 0 ) THEN 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 ! Test if "no night" point 755 IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 756 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 757 & / icount_sstnight(ji,jj) 758 ELSE 759 sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 760 ENDIF 761 END DO 762 END DO 763 ENDIF 764 765 ENDIF 766 767 ! Get the data for interpolation 768 769 ALLOCATE( & 770 & igrdi(2,2,isst), & 771 & igrdj(2,2,isst), & 772 & zglam(2,2,isst), & 773 & zgphi(2,2,isst), & 774 & zmask(2,2,isst), & 775 & zsstl(2,2,isst) & 776 & ) 777 778 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 779 iobs = jobs - sstdatqc%nsurfup 780 igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 781 igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 782 igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 783 igrdj(1,2,iobs) = sstdatqc%mj(jobs) 784 igrdi(2,1,iobs) = sstdatqc%mi(jobs) 785 igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 786 igrdi(2,2,iobs) = sstdatqc%mi(jobs) 787 igrdj(2,2,iobs) = sstdatqc%mj(jobs) 788 END DO 789 790 CALL obs_int_comm_2d( 2, 2, isst, & 791 & igrdi, igrdj, glamt, zglam ) 792 CALL obs_int_comm_2d( 2, 2, isst, & 793 & igrdi, igrdj, gphit, zgphi ) 794 CALL obs_int_comm_2d( 2, 2, isst, & 795 & igrdi, igrdj, psstmask, zmask ) 796 CALL obs_int_comm_2d( 2, 2, isst, & 797 & igrdi, igrdj, psstn, zsstl ) 798 799 ! At the end of the day get interpolated means 800 IF ( idayend == 0 .AND. ld_nightav ) THEN 801 802 ALLOCATE( & 803 & zsstm(2,2,isst) & 804 & ) 805 806 CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 807 & sstdatqc%vdmean(:,:), zsstm ) 808 809 ENDIF 810 811 ! Loop over observations 812 813 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 814 815 iobs = jobs - sstdatqc%nsurfup 816 817 IF ( kt /= sstdatqc%mstp(jobs) ) THEN 818 819 IF(lwp) THEN 820 WRITE(numout,*) 821 WRITE(numout,*) ' E R R O R : Observation', & 822 & ' time step is not consistent with the', & 823 & ' model time step' 824 WRITE(numout,*) ' =========' 825 WRITE(numout,*) 826 WRITE(numout,*) ' Record = ', jobs, & 827 & ' kt = ', kt, & 828 & ' mstp = ', sstdatqc%mstp(jobs), & 829 & ' ntyp = ', sstdatqc%ntyp(jobs) 830 ENDIF 831 CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 832 833 ENDIF 834 835 zlam = sstdatqc%rlam(jobs) 836 zphi = sstdatqc%rphi(jobs) 837 838 ! Get weights to interpolate the model SST to the observation point 839 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 840 & zglam(:,:,iobs), zgphi(:,:,iobs), & 841 & zmask(:,:,iobs), zweig, zobsmask ) 842 843 ! Interpolate the model SST to the observation point 844 845 IF ( ld_nightav ) THEN 846 847 IF ( idayend == 0 ) THEN 848 ! Daily averaged/diurnal cycle of SST data 849 CALL obs_int_h2d( 1, 1, & 850 & zweig, zsstm(:,:,iobs), zext ) 851 ELSE 852 CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // & 853 & ' number of night SST data should' // & 854 & ' only occur at the end of a given day' ) 855 ENDIF 856 857 ELSE 858 859 CALL obs_int_h2d( 1, 1, & 860 & zweig, zsstl(:,:,iobs), zext ) 861 862 ENDIF 863 sstdatqc%rmod(jobs,1) = zext(1) 864 865 END DO 866 867 ! Deallocate the data for interpolation 868 DEALLOCATE( & 869 & igrdi, & 870 & igrdj, & 871 & zglam, & 872 & zgphi, & 873 & zmask, & 874 & zsstl & 875 & ) 876 877 ! At the end of the day also get interpolated means 878 IF ( idayend == 0 .AND. ld_nightav ) THEN 1335 ! At the end of the day also deallocate night-time mean array 1336 IF ( idayend == 0 .AND. ldnightav ) THEN 879 1337 DEALLOCATE( & 880 & zs stm &1338 & zsurfm & 881 1339 & ) 882 1340 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 1341 1342 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1343 1344 END SUBROUTINE obs_surf_opt 1440 1345 1441 1346 END MODULE obs_oper 1442
Note: See TracChangeset
for help on using the changeset viewer.
