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