Changeset 10120 for branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
- Timestamp:
- 2018-09-12T19:12:15+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r6486 r10120 7 7 8 8 !!---------------------------------------------------------------------- 9 !! obs_pro_opt : Compute the model counterpart of temperature and 10 !! salinity observations from profiles 11 !! obs_sla_opt : Compute the model counterpart of sea level anomaly 12 !! observations 13 !! obs_sst_opt : Compute the model counterpart of sea surface temperature 14 !! observations 15 !! obs_sss_opt : Compute the model counterpart of sea surface salinity 16 !! observations 17 !! obs_seaice_opt : Compute the model counterpart of sea ice concentration 18 !! observations 19 !! 20 !! obs_vel_opt : Compute the model counterpart of zonal and meridional 21 !! components of velocity from observations. 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 !! obs_surf_opt : Compute the model counterpart of surface data 22 11 !!---------------------------------------------------------------------- 23 12 24 !! * Modules used 13 !! * Modules used 25 14 USE par_kind, ONLY : & ! Precision variables 26 15 & wp 27 16 USE in_out_manager ! I/O manager 28 17 USE obs_inter_sup ! Interpolation support 29 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs ervationpt18 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 30 19 & obs_int_h2d, & 31 20 & obs_int_h2d_init 32 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 33 26 & obs_int_z1d, & 34 27 & obs_int_z1d_spl 35 USE obs_const, ONLY : &36 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 37 30 USE dom_oce, ONLY : & 38 & glamt, glam u, glamv, &39 & gphit, gphi u, gphiv40 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 41 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 37 USE obs_grid, ONLY : & 38 & obs_level_search 42 39 43 40 IMPLICIT NONE … … 46 43 PRIVATE 47 44 48 PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations 49 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 50 & obs_sst_opt, & ! Compute the model counterpart of SST observations 51 & obs_sss_opt, & ! Compute the model counterpart of SSS observations 52 & obs_seaice_opt, & 53 & obs_vel_opt ! Compute the model counterpart of velocity profile data 54 55 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_surf_opt ! Compute the model counterpart of surface obs 47 48 INTEGER, PARAMETER, PUBLIC :: & 49 & imaxavtypes = 20 ! Max number of daily avgd obs types 56 50 57 51 !!---------------------------------------------------------------------- … … 61 55 !!---------------------------------------------------------------------- 62 56 57 !! * Substitutions 58 # include "domzgr_substitute.h90" 63 59 CONTAINS 64 60 65 SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 66 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 67 & kdailyavtypes ) 61 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 63 & kit000, kdaystp, & 64 & pvar1, pvar2, pgdept, pgdepw, & 65 & pmask1, pmask2, & 66 & plam1, plam2, pphi1, pphi2, & 67 & k1dint, k2dint, kdailyavtypes ) 68 68 69 !!----------------------------------------------------------------------- 69 70 !! … … 78 79 !! 79 80 !! First, a vertical profile of horizontally interpolated model 80 !! now temperatures is computed at the obs (lon, lat) point.81 !! now values is computed at the obs (lon, lat) point. 81 82 !! Several horizontal interpolation schemes are available: 82 83 !! - distance-weighted (great circle) (k2dint = 0) … … 86 87 !! - polynomial (quadrilateral grid) (k2dint = 4) 87 88 !! 88 !! Next, the vertical temperatureprofile is interpolated to the89 !! Next, the vertical profile is interpolated to the 89 90 !! data depth points. Two vertical interpolation schemes are 90 91 !! available: … … 96 97 !! routine. 97 98 !! 98 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is99 !! If the logical is switched on, the model equivalent is 99 100 !! a daily mean model temperature field. So, we first compute 100 101 !! the mean, then interpolate only at the end of the day. 101 102 !! 102 !! Note: thein situ temperature observations must be converted103 !! Note: in situ temperature observations must be converted 103 104 !! to potential temperature (the model variable) prior to 104 105 !! assimilation. 105 !!??????????????????????????????????????????????????????????????106 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???107 !!??????????????????????????????????????????????????????????????108 106 !! 109 107 !! ** Action : … … 115 113 !! ! 07-01 (K. Mogensen) Merge of temperature and salinity 116 114 !! ! 07-03 (K. Mogensen) General handling of profiles 115 !! ! 15-02 (M. Martin) Combined routine for all profile types 116 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 117 117 !!----------------------------------------------------------------------- 118 118 119 119 !! * Modules used 120 120 USE obs_profiles_def ! Definition of storage space for profile obs. … … 123 123 124 124 !! * Arguments 125 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 126 INTEGER, INTENT(IN) :: kt ! Time step 127 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 125 TYPE(obs_prof), INTENT(INOUT) :: & 126 & prodatqc ! Subset of profile data passing QC 127 INTEGER, INTENT(IN) :: kt ! Time step 128 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 128 129 INTEGER, INTENT(IN) :: kpj 129 130 INTEGER, INTENT(IN) :: kpk 130 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step131 132 INTEGER, INTENT(IN) :: k1dint 133 INTEGER, INTENT(IN) :: k2dint 134 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day131 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 132 ! (kit000-1 = restart time) 133 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 134 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 135 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 135 136 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 136 & ptn, & ! Model temperature field 137 & psn, & ! Model salinity field 138 & ptmask ! Land-sea mask 139 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 140 & pgdept ! Model array of depth levels 137 & pvar1, & ! Model field 1 138 & pvar2, & ! Model field 2 139 & pmask1, & ! Land-sea mask 1 140 & pmask2 ! Land-sea mask 2 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 142 & plam1, & ! Model longitudes for variable 1 143 & plam2, & ! Model longitudes for variable 2 144 & pphi1, & ! Model latitudes for variable 1 145 & pphi2 ! Model latitudes for variable 2 146 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 147 & pgdept, & ! Model array of depth T levels 148 & pgdepw ! Model array of depth W levels 141 149 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 142 & kdailyavtypes! Types for daily averages 150 & kdailyavtypes ! Types for daily averages 151 143 152 !! * Local declarations 144 153 INTEGER :: ji … … 152 161 INTEGER :: iend 153 162 INTEGER :: iobs 163 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 164 INTEGER :: inum_obs 154 165 INTEGER, DIMENSION(imaxavtypes) :: & 155 166 & idailyavtypes 167 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 168 & igrdi1, & 169 & igrdi2, & 170 & igrdj1, & 171 & igrdj2 172 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 173 156 174 REAL(KIND=wp) :: zlam 157 175 REAL(KIND=wp) :: zphi 158 176 REAL(KIND=wp) :: zdaystp 159 177 REAL(KIND=wp), DIMENSION(kpk) :: & 160 & zobsmask, & 178 & zobsmask1, & 179 & zobsmask2, & 161 180 & zobsk, & 162 181 & zobs2k 163 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 182 REAL(KIND=wp), DIMENSION(2,2,1) :: & 183 & zweig1, & 184 & zweig2, & 164 185 & zweig 165 186 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 166 & zmask, & 167 & zintt, & 168 & zints, & 169 & zinmt, & 170 & zinms 187 & zmask1, & 188 & zmask2, & 189 & zint1, & 190 & zint2, & 191 & zinm1, & 192 & zinm2, & 193 & zgdept, & 194 & zgdepw 171 195 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 172 & zglam, & 173 & zgphi 174 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 175 & igrdi, & 176 & igrdj 196 & zglam1, & 197 & zglam2, & 198 & zgphi1, & 199 & zgphi2 200 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 201 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 202 203 LOGICAL :: ld_dailyav 177 204 178 205 !------------------------------------------------------------------------ 179 206 ! Local initialization 180 207 !------------------------------------------------------------------------ 181 ! ...Record and data counters208 ! Record and data counters 182 209 inrc = kt - kit000 + 2 183 210 ipro = prodatqc%npstp(inrc) 184 211 185 212 ! Daily average types 213 ld_dailyav = .FALSE. 186 214 IF ( PRESENT(kdailyavtypes) ) THEN 187 215 idailyavtypes(:) = kdailyavtypes(:) 216 IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 188 217 ELSE 189 218 idailyavtypes(:) = -1 190 219 ENDIF 191 220 192 ! Initialize daily mean for first timestep 221 ! Daily means are calculated for values over timesteps: 222 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 193 223 idayend = MOD( kt - kit000 + 1, kdaystp ) 194 224 195 ! Added kt == 0 test to catch restart case 196 IF ( idayend == 1 .OR. kt == 0) THEN 197 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 225 IF ( ld_dailyav ) THEN 226 227 ! Initialize daily mean for first timestep of the day 228 IF ( idayend == 1 .OR. kt == 0 ) THEN 229 DO jk = 1, jpk 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 prodatqc%vdmean(ji,jj,jk,1) = 0.0 233 prodatqc%vdmean(ji,jj,jk,2) = 0.0 234 END DO 235 END DO 236 END DO 237 ENDIF 238 198 239 DO jk = 1, jpk 199 240 DO jj = 1, jpj 200 241 DO ji = 1, jpi 201 prodatqc%vdmean(ji,jj,jk,1) = 0.0 202 prodatqc%vdmean(ji,jj,jk,2) = 0.0 242 ! Increment field 1 for computing daily mean 243 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 244 & + pvar1(ji,jj,jk) 245 ! Increment field 2 for computing daily mean 246 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 247 & + pvar2(ji,jj,jk) 203 248 END DO 204 249 END DO 205 250 END DO 206 ENDIF 207 208 DO jk = 1, jpk 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 ! Increment the temperature field for computing daily mean 212 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 213 & + ptn(ji,jj,jk) 214 ! Increment the salinity field for computing daily mean 215 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 216 & + psn(ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 ! Compute the daily mean at the end of day 222 zdaystp = 1.0 / REAL( kdaystp ) 223 IF ( idayend == 0 ) THEN 224 DO jk = 1, jpk 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 228 & * zdaystp 229 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 230 & * zdaystp 251 252 ! Compute the daily mean at the end of day 253 zdaystp = 1.0 / REAL( kdaystp ) 254 IF ( idayend == 0 ) THEN 255 IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 256 CALL FLUSH(numout) 257 DO jk = 1, jpk 258 DO jj = 1, jpj 259 DO ji = 1, jpi 260 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 261 & * zdaystp 262 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 263 & * zdaystp 264 END DO 231 265 END DO 232 266 END DO 233 END DO 267 ENDIF 268 234 269 ENDIF 235 270 236 271 ! Get the data for interpolation 237 272 ALLOCATE( & 238 & igrdi(2,2,ipro), & 239 & igrdj(2,2,ipro), & 240 & zglam(2,2,ipro), & 241 & zgphi(2,2,ipro), & 242 & zmask(2,2,kpk,ipro), & 243 & zintt(2,2,kpk,ipro), & 244 & zints(2,2,kpk,ipro) & 273 & igrdi1(2,2,ipro), & 274 & igrdi2(2,2,ipro), & 275 & igrdj1(2,2,ipro), & 276 & igrdj2(2,2,ipro), & 277 & zglam1(2,2,ipro), & 278 & zglam2(2,2,ipro), & 279 & zgphi1(2,2,ipro), & 280 & zgphi2(2,2,ipro), & 281 & zmask1(2,2,kpk,ipro), & 282 & zmask2(2,2,kpk,ipro), & 283 & zint1(2,2,kpk,ipro), & 284 & zint2(2,2,kpk,ipro), & 285 & zgdept(2,2,kpk,ipro), & 286 & zgdepw(2,2,kpk,ipro) & 245 287 & ) 246 288 247 289 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 248 290 iobs = jobs - prodatqc%nprofup 249 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 250 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 251 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 252 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 253 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 254 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 255 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 256 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 291 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 292 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 293 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 294 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 295 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 296 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 297 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 298 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 299 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 300 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 301 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 302 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 303 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 304 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 305 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 306 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 257 307 END DO 258 308 259 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 260 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 261 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 262 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 263 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 309 ! Initialise depth arrays 310 zgdept(:,:,:,:) = 0.0 311 zgdepw(:,:,:,:) = 0.0 312 313 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 314 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 315 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 316 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 317 318 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 319 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 320 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 322 323 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 324 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 264 325 265 326 ! At the end of the day also get interpolated means 266 IF ( idayend == 0 ) THEN327 IF ( ld_dailyav .AND. idayend == 0 ) THEN 267 328 268 329 ALLOCATE( & 269 & zinm t(2,2,kpk,ipro), &270 & zinm s(2,2,kpk,ipro) &330 & zinm1(2,2,kpk,ipro), & 331 & zinm2(2,2,kpk,ipro) & 271 332 & ) 272 333 273 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &274 & prodatqc%vdmean(:,:,:,1), zinm t)275 CALL obs_int_comm_3d( 2, 2, ipro, kp k, igrdi, igrdj, &276 & prodatqc%vdmean(:,:,:,2), zinm s)334 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 335 & prodatqc%vdmean(:,:,:,1), zinm1 ) 336 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 337 & prodatqc%vdmean(:,:,:,2), zinm2 ) 277 338 278 339 ENDIF 279 340 341 ! Return if no observations to process 342 ! Has to be done after comm commands to ensure processors 343 ! stay in sync 344 IF ( ipro == 0 ) RETURN 345 280 346 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 281 347 … … 283 349 284 350 IF ( kt /= prodatqc%mstp(jobs) ) THEN 285 351 286 352 IF(lwp) THEN 287 353 WRITE(numout,*) … … 298 364 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 299 365 ENDIF 300 366 301 367 zlam = prodatqc%rlam(jobs) 302 368 zphi = prodatqc%rphi(jobs) 369 370 ! Horizontal weights 371 ! Masked values are calculated later. 372 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 373 374 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 375 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 376 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 377 378 ENDIF 379 380 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 381 382 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 383 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 384 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 385 386 ENDIF 387 388 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 389 390 zobsk(:) = obfillflt 391 392 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 393 394 IF ( idayend == 0 ) THEN 395 ! Daily averaged data 396 397 ! vertically interpolate all 4 corners 398 ista = prodatqc%npvsta(jobs,1) 399 iend = prodatqc%npvend(jobs,1) 400 inum_obs = iend - ista + 1 401 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 402 403 DO iin=1,2 404 DO ijn=1,2 405 406 IF ( k1dint == 1 ) THEN 407 CALL obs_int_z1d_spl( kpk, & 408 & zinm1(iin,ijn,:,iobs), & 409 & zobs2k, zgdept(iin,ijn,:,iobs), & 410 & zmask1(iin,ijn,:,iobs)) 411 ENDIF 412 413 CALL obs_level_search(kpk, & 414 & zgdept(iin,ijn,:,iobs), & 415 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 416 & iv_indic) 417 418 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 419 & prodatqc%var(1)%vdep(ista:iend), & 420 & zinm1(iin,ijn,:,iobs), & 421 & zobs2k, interp_corner(iin,ijn,:), & 422 & zgdept(iin,ijn,:,iobs), & 423 & zmask1(iin,ijn,:,iobs)) 424 425 ENDDO 426 ENDDO 427 428 ENDIF !idayend 429 430 ELSE 431 432 ! Point data 433 434 ! vertically interpolate all 4 corners 435 ista = prodatqc%npvsta(jobs,1) 436 iend = prodatqc%npvend(jobs,1) 437 inum_obs = iend - ista + 1 438 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 439 DO iin=1,2 440 DO ijn=1,2 441 442 IF ( k1dint == 1 ) THEN 443 CALL obs_int_z1d_spl( kpk, & 444 & zint1(iin,ijn,:,iobs),& 445 & zobs2k, zgdept(iin,ijn,:,iobs), & 446 & zmask1(iin,ijn,:,iobs)) 447 448 ENDIF 449 450 CALL obs_level_search(kpk, & 451 & zgdept(iin,ijn,:,iobs),& 452 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 453 & iv_indic) 454 455 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 456 & prodatqc%var(1)%vdep(ista:iend), & 457 & zint1(iin,ijn,:,iobs), & 458 & zobs2k,interp_corner(iin,ijn,:), & 459 & zgdept(iin,ijn,:,iobs), & 460 & zmask1(iin,ijn,:,iobs) ) 303 461 304 ! Horizontal weights and vertical mask 305 306 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 307 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 308 309 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 310 & zglam(:,:,iobs), zgphi(:,:,iobs), & 311 & zmask(:,:,:,iobs), zweig, zobsmask ) 312 313 ENDIF 314 315 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 462 ENDDO 463 ENDDO 464 465 ENDIF 466 467 !------------------------------------------------------------- 468 ! Compute the horizontal interpolation for every profile level 469 !------------------------------------------------------------- 470 471 DO ikn=1,inum_obs 472 iend=ista+ikn-1 473 474 zweig(:,:,1) = 0._wp 475 476 ! This code forces the horizontal weights to be 477 ! zero IF the observation is below the bottom of the 478 ! corners of the interpolation nodes, Or if it is in 479 ! the mask. This is important for observations near 480 ! steep bathymetry 481 DO iin=1,2 482 DO ijn=1,2 483 484 depth_loop1: DO ik=kpk,2,-1 485 IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN 486 487 zweig(iin,ijn,1) = & 488 & zweig1(iin,ijn,1) * & 489 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 490 & - prodatqc%var(1)%vdep(iend)),0._wp) 491 492 EXIT depth_loop1 493 494 ENDIF 495 496 ENDDO depth_loop1 497 498 ENDDO 499 ENDDO 500 501 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 502 & prodatqc%var(1)%vmod(iend:iend) ) 503 504 ! Set QC flag for any observations found below the bottom 505 ! needed as the check here is more strict than that in obs_prep 506 IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 507 508 ENDDO 509 510 DEALLOCATE(interp_corner,iv_indic) 511 512 ENDIF 513 514 ! For the second variable 515 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 316 516 317 517 zobsk(:) = obfillflt 318 518 319 519 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 320 520 321 521 IF ( idayend == 0 ) THEN 522 ! Daily averaged data 523 524 ! vertically interpolate all 4 corners 525 ista = prodatqc%npvsta(jobs,2) 526 iend = prodatqc%npvend(jobs,2) 527 inum_obs = iend - ista + 1 528 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 529 530 DO iin=1,2 531 DO ijn=1,2 532 533 IF ( k1dint == 1 ) THEN 534 CALL obs_int_z1d_spl( kpk, & 535 & zinm2(iin,ijn,:,iobs), & 536 & zobs2k, zgdept(iin,ijn,:,iobs), & 537 & zmask2(iin,ijn,:,iobs)) 538 ENDIF 539 540 CALL obs_level_search(kpk, & 541 & zgdept(iin,ijn,:,iobs), & 542 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 543 & iv_indic) 544 545 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 546 & prodatqc%var(2)%vdep(ista:iend), & 547 & zinm2(iin,ijn,:,iobs), & 548 & zobs2k, interp_corner(iin,ijn,:), & 549 & zgdept(iin,ijn,:,iobs), & 550 & zmask2(iin,ijn,:,iobs)) 551 552 ENDDO 553 ENDDO 554 555 ENDIF !idayend 556 557 ELSE 558 559 ! Point data 560 561 ! vertically interpolate all 4 corners 562 ista = prodatqc%npvsta(jobs,2) 563 iend = prodatqc%npvend(jobs,2) 564 inum_obs = iend - ista + 1 565 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 566 DO iin=1,2 567 DO ijn=1,2 568 569 IF ( k1dint == 1 ) THEN 570 CALL obs_int_z1d_spl( kpk, & 571 & zint2(iin,ijn,:,iobs),& 572 & zobs2k, zgdept(iin,ijn,:,iobs), & 573 & zmask2(iin,ijn,:,iobs)) 574 575 ENDIF 576 577 CALL obs_level_search(kpk, & 578 & zgdept(iin,ijn,:,iobs),& 579 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 580 & iv_indic) 581 582 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 583 & prodatqc%var(2)%vdep(ista:iend), & 584 & zint2(iin,ijn,:,iobs), & 585 & zobs2k,interp_corner(iin,ijn,:), & 586 & zgdept(iin,ijn,:,iobs), & 587 & zmask2(iin,ijn,:,iobs) ) 588 589 ENDDO 590 ENDDO 591 592 ENDIF 593 594 !------------------------------------------------------------- 595 ! Compute the horizontal interpolation for every profile level 596 !------------------------------------------------------------- 597 598 DO ikn=1,inum_obs 599 iend=ista+ikn-1 322 600 323 ! Daily averaged moored buoy (MRB) data 324 325 CALL obs_int_h2d( kpk, kpk, & 326 & zweig, zinmt(:,:,:,iobs), zobsk ) 327 328 329 ELSE 330 331 CALL ctl_stop( ' A nonzero' // & 332 & ' number of profile T BUOY data should' // & 333 & ' only occur at the end of a given day' ) 334 335 ENDIF 336 337 ELSE 338 339 ! Point data 340 341 CALL obs_int_h2d( kpk, kpk, & 342 & zweig, zintt(:,:,:,iobs), zobsk ) 343 344 ENDIF 345 346 !------------------------------------------------------------- 347 ! Compute vertical second-derivative of the interpolating 348 ! polynomial at obs points 349 !------------------------------------------------------------- 350 351 IF ( k1dint == 1 ) THEN 352 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 353 & pgdept, zobsmask ) 354 ENDIF 355 356 !----------------------------------------------------------------- 357 ! Vertical interpolation to the observation point 358 !----------------------------------------------------------------- 359 ista = prodatqc%npvsta(jobs,1) 360 iend = prodatqc%npvend(jobs,1) 361 CALL obs_int_z1d( kpk, & 362 & prodatqc%var(1)%mvk(ista:iend), & 363 & k1dint, iend - ista + 1, & 364 & prodatqc%var(1)%vdep(ista:iend), & 365 & zobsk, zobs2k, & 366 & prodatqc%var(1)%vmod(ista:iend), & 367 & pgdept, zobsmask ) 368 369 ENDIF 370 371 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 372 373 zobsk(:) = obfillflt 374 375 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 376 377 IF ( idayend == 0 ) THEN 378 379 ! Daily averaged moored buoy (MRB) data 380 381 CALL obs_int_h2d( kpk, kpk, & 382 & zweig, zinms(:,:,:,iobs), zobsk ) 383 384 ELSE 385 386 CALL ctl_stop( ' A nonzero' // & 387 & ' number of profile S BUOY data should' // & 388 & ' only occur at the end of a given day' ) 389 390 ENDIF 391 392 ELSE 393 394 ! Point data 395 396 CALL obs_int_h2d( kpk, kpk, & 397 & zweig, zints(:,:,:,iobs), zobsk ) 398 399 ENDIF 400 401 402 !------------------------------------------------------------- 403 ! Compute vertical second-derivative of the interpolating 404 ! polynomial at obs points 405 !------------------------------------------------------------- 406 407 IF ( k1dint == 1 ) THEN 408 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 409 & pgdept, zobsmask ) 410 ENDIF 411 412 !---------------------------------------------------------------- 413 ! Vertical interpolation to the observation point 414 !---------------------------------------------------------------- 415 ista = prodatqc%npvsta(jobs,2) 416 iend = prodatqc%npvend(jobs,2) 417 CALL obs_int_z1d( kpk, & 418 & prodatqc%var(2)%mvk(ista:iend),& 419 & k1dint, iend - ista + 1, & 420 & prodatqc%var(2)%vdep(ista:iend),& 421 & zobsk, zobs2k, & 422 & prodatqc%var(2)%vmod(ista:iend),& 423 & pgdept, zobsmask ) 424 425 ENDIF 426 427 END DO 601 zweig(:,:,1) = 0._wp 602 603 ! This code forces the horizontal weights to be 604 ! zero IF the observation is below the bottom of the 605 ! corners of the interpolation nodes, Or if it is in 606 ! the mask. This is important for observations near 607 ! steep bathymetry 608 DO iin=1,2 609 DO ijn=1,2 610 611 depth_loop2: DO ik=kpk,2,-1 612 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 613 614 zweig(iin,ijn,1) = & 615 & zweig2(iin,ijn,1) * & 616 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 617 & - prodatqc%var(2)%vdep(iend)),0._wp) 618 619 EXIT depth_loop2 620 621 ENDIF 622 623 ENDDO depth_loop2 624 625 ENDDO 626 ENDDO 627 628 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 629 & prodatqc%var(2)%vmod(iend:iend) ) 630 631 ! Set QC flag for any observations found below the bottom 632 ! needed as the check here is more strict than that in obs_prep 633 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 428 634 635 ENDDO 636 637 DEALLOCATE(interp_corner,iv_indic) 638 639 ENDIF 640 641 ENDDO 642 429 643 ! Deallocate the data for interpolation 430 644 DEALLOCATE( & 431 & igrdi, & 432 & igrdj, & 433 & zglam, & 434 & zgphi, & 435 & zmask, & 436 & zintt, & 437 & zints & 645 & igrdi1, & 646 & igrdi2, & 647 & igrdj1, & 648 & igrdj2, & 649 & zglam1, & 650 & zglam2, & 651 & zgphi1, & 652 & zgphi2, & 653 & zmask1, & 654 & zmask2, & 655 & zint1, & 656 & zint2, & 657 & zgdept, & 658 & zgdepw & 438 659 & ) 660 439 661 ! At the end of the day also get interpolated means 440 IF ( idayend == 0 ) THEN662 IF ( ld_dailyav .AND. idayend == 0 ) THEN 441 663 DEALLOCATE( & 442 & zinm t, &443 & zinm s&664 & zinm1, & 665 & zinm2 & 444 666 & ) 445 667 ENDIF 446 668 447 669 prodatqc%nprofup = prodatqc%nprofup + ipro 448 449 END SUBROUTINE obs_pro_opt 450 451 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 & psshn, psshmask, k2dint ) 670 671 END SUBROUTINE obs_prof_opt 672 673 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 674 & kit000, kdaystp, psurf, psurfmask, & 675 & k2dint, ldnightav, plamscl, pphiscl, & 676 & lindegrees ) 677 453 678 !!----------------------------------------------------------------------- 454 679 !! 455 !! *** ROUTINE obs_s la_opt ***456 !! 457 !! ** Purpose : Compute the model counterpart of s ea level anomaly680 !! *** ROUTINE obs_surf_opt *** 681 !! 682 !! ** Purpose : Compute the model counterpart of surface 458 683 !! data by interpolating from the model grid to the 459 684 !! observation point. … … 462 687 !! the model values at the corners of the surrounding grid box. 463 688 !! 464 !! The n ow model SSHis first computed at the obs (lon, lat) point.689 !! The new model value is first computed at the obs (lon, lat) point. 465 690 !! 466 691 !! Several horizontal interpolation schemes are available: … … 470 695 !! - bilinear (quadrilateral grid) (k2dint = 3) 471 696 !! - polynomial (quadrilateral grid) (k2dint = 4) 472 !! 473 !! The sea level anomaly at the observation points is then computed 474 !! by removing a mean dynamic topography (defined at the obs. point). 697 !! 698 !! Two horizontal averaging schemes are also available: 699 !! - weighted radial footprint (k2dint = 5) 700 !! - weighted rectangular footprint (k2dint = 6) 701 !! 475 702 !! 476 703 !! ** Action : … … 478 705 !! History : 479 706 !! ! 07-03 (A. Weaver) 707 !! ! 15-02 (M. Martin) Combined routine for surface types 708 !! ! 17-03 (M. Martin) Added horizontal averaging options 480 709 !!----------------------------------------------------------------------- 481 710 482 711 !! * Modules used 483 712 USE obs_surf_def ! Definition of storage space for surface observations … … 486 715 487 716 !! * Arguments 488 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening 489 INTEGER, INTENT(IN) :: kt ! Time step 490 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 717 TYPE(obs_surf), INTENT(INOUT) :: & 718 & surfdataqc ! Subset of surface data passing QC 719 INTEGER, INTENT(IN) :: kt ! Time step 720 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 491 721 INTEGER, INTENT(IN) :: kpj 492 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 493 ! (kit000-1 = restart time) 494 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 495 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 496 & psshn, & ! Model SSH field 497 & psshmask ! Land-sea mask 498 722 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 723 ! (kit000-1 = restart time) 724 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 725 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 726 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 727 & psurf, & ! Model surface field 728 & psurfmask ! Land-sea mask 729 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 730 REAL(KIND=wp), INTENT(IN) :: & 731 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 732 & pphiscl ! This is the full width (rather than half-width) 733 LOGICAL, INTENT(IN) :: & 734 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 735 499 736 !! * Local declarations 500 737 INTEGER :: ji … … 502 739 INTEGER :: jobs 503 740 INTEGER :: inrc 504 INTEGER :: is la741 INTEGER :: isurf 505 742 INTEGER :: iobs 506 REAL(KIND=wp) :: zlam 507 REAL(KIND=wp) :: zphi 508 REAL(KIND=wp) :: zext(1), zobsmask(1) 509 REAL(kind=wp), DIMENSION(2,2,1) :: & 510 & zweig 743 INTEGER :: imaxifp, imaxjfp 744 INTEGER :: imodi, imodj 745 INTEGER :: idayend 746 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 747 & igrdi, & 748 & igrdj, & 749 & igrdip1, & 750 & igrdjp1 751 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 752 & icount_night, & 753 & imask_night 754 REAL(wp) :: zlam 755 REAL(wp) :: zphi 756 REAL(wp), DIMENSION(1) :: zext, zobsmask 757 REAL(wp) :: zdaystp 511 758 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 512 & zmask, & 513 & zsshl, & 514 & zglam, & 515 & zgphi 516 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 517 & igrdi, & 518 & igrdj 759 & zweig, & 760 & zmask, & 761 & zsurf, & 762 & zsurfm, & 763 & zsurftmp, & 764 & zglam, & 765 & zgphi, & 766 & zglamf, & 767 & zgphif 768 769 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 770 & zintmp, & 771 & zouttmp, & 772 & zmeanday ! to compute model sst in region of 24h daylight (pole) 519 773 520 774 !------------------------------------------------------------------------ 521 775 ! Local initialization 522 776 !------------------------------------------------------------------------ 523 ! ...Record and data counters777 ! Record and data counters 524 778 inrc = kt - kit000 + 2 525 isla = sladatqc%nsstp(inrc) 779 isurf = surfdataqc%nsstp(inrc) 780 781 ! Work out the maximum footprint size for the 782 ! interpolation/averaging in model grid-points - has to be even. 783 784 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 785 786 787 IF ( ldnightav ) THEN 788 789 ! Initialize array for night mean 790 IF ( kt == 0 ) THEN 791 ALLOCATE ( icount_night(kpi,kpj) ) 792 ALLOCATE ( imask_night(kpi,kpj) ) 793 ALLOCATE ( zintmp(kpi,kpj) ) 794 ALLOCATE ( zouttmp(kpi,kpj) ) 795 ALLOCATE ( zmeanday(kpi,kpj) ) 796 nday_qsr = -1 ! initialisation flag for nbc_dcy 797 ENDIF 798 799 ! Night-time means are calculated for night-time values over timesteps: 800 ! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 801 idayend = MOD( kt - kit000 + 1, kdaystp ) 802 803 ! Initialize night-time mean for first timestep of the day 804 IF ( idayend == 1 .OR. kt == 0 ) THEN 805 DO jj = 1, jpj 806 DO ji = 1, jpi 807 surfdataqc%vdmean(ji,jj) = 0.0 808 zmeanday(ji,jj) = 0.0 809 icount_night(ji,jj) = 0 810 END DO 811 END DO 812 ENDIF 813 814 zintmp(:,:) = 0.0 815 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 816 imask_night(:,:) = INT( zouttmp(:,:) ) 817 818 DO jj = 1, jpj 819 DO ji = 1, jpi 820 ! Increment the temperature field for computing night mean and counter 821 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 822 & + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 823 zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj) 824 icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 825 END DO 826 END DO 827 828 ! Compute the night-time mean at the end of the day 829 zdaystp = 1.0 / REAL( kdaystp ) 830 IF ( idayend == 0 ) THEN 831 IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 832 DO jj = 1, jpj 833 DO ji = 1, jpi 834 ! Test if "no night" point 835 IF ( icount_night(ji,jj) > 0 ) THEN 836 surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 837 & / REAL( icount_night(ji,jj) ) 838 ELSE 839 !At locations where there is no night (e.g. poles), 840 ! calculate daily mean instead of night-time mean. 841 surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 842 ENDIF 843 END DO 844 END DO 845 ENDIF 846 847 ENDIF 526 848 527 849 ! Get the data for interpolation 528 850 529 851 ALLOCATE( & 530 & igrdi(2,2,isla), & 531 & igrdj(2,2,isla), & 532 & zglam(2,2,isla), & 533 & zgphi(2,2,isla), & 534 & zmask(2,2,isla), & 535 & zsshl(2,2,isla) & 852 & zweig(imaxifp,imaxjfp,1), & 853 & igrdi(imaxifp,imaxjfp,isurf), & 854 & igrdj(imaxifp,imaxjfp,isurf), & 855 & zglam(imaxifp,imaxjfp,isurf), & 856 & zgphi(imaxifp,imaxjfp,isurf), & 857 & zmask(imaxifp,imaxjfp,isurf), & 858 & zsurf(imaxifp,imaxjfp,isurf), & 859 & zsurftmp(imaxifp,imaxjfp,isurf), & 860 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 861 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 862 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 863 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 536 864 & ) 537 538 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 539 iobs = jobs - sladatqc%nsurfup 540 igrdi(1,1,iobs) = sladatqc%mi(jobs)-1 541 igrdj(1,1,iobs) = sladatqc%mj(jobs)-1 542 igrdi(1,2,iobs) = sladatqc%mi(jobs)-1 543 igrdj(1,2,iobs) = sladatqc%mj(jobs) 544 igrdi(2,1,iobs) = sladatqc%mi(jobs) 545 igrdj(2,1,iobs) = sladatqc%mj(jobs)-1 546 igrdi(2,2,iobs) = sladatqc%mi(jobs) 547 igrdj(2,2,iobs) = sladatqc%mj(jobs) 865 866 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 867 iobs = jobs - surfdataqc%nsurfup 868 DO ji = 0, imaxifp 869 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 870 871 !Deal with wrap around in longitude 872 IF ( imodi < 1 ) imodi = imodi + jpiglo 873 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 874 875 DO jj = 0, imaxjfp 876 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 877 !If model values are out of the domain to the north/south then 878 !set them to be the edge of the domain 879 IF ( imodj < 1 ) imodj = 1 880 IF ( imodj > jpjglo ) imodj = jpjglo 881 882 igrdip1(ji+1,jj+1,iobs) = imodi 883 igrdjp1(ji+1,jj+1,iobs) = imodj 884 885 IF ( ji >= 1 .AND. jj >= 1 ) THEN 886 igrdi(ji,jj,iobs) = imodi 887 igrdj(ji,jj,iobs) = imodj 888 ENDIF 889 890 END DO 891 END DO 548 892 END DO 549 893 550 CALL obs_int_comm_2d( 2, 2, isla, &894 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 551 895 & igrdi, igrdj, glamt, zglam ) 552 CALL obs_int_comm_2d( 2, 2, isla, &896 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 553 897 & igrdi, igrdj, gphit, zgphi ) 554 CALL obs_int_comm_2d( 2, 2, isla, & 555 & igrdi, igrdj, psshmask, zmask ) 556 CALL obs_int_comm_2d( 2, 2, isla, & 557 & igrdi, igrdj, psshn, zsshl ) 898 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 899 & igrdi, igrdj, psurfmask, zmask ) 900 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 901 & igrdi, igrdj, psurf, zsurf ) 902 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 903 & igrdip1, igrdjp1, glamf, zglamf ) 904 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 905 & igrdip1, igrdjp1, gphif, zgphif ) 906 907 ! At the end of the day get interpolated means 908 IF ( idayend == 0 .AND. ldnightav ) THEN 909 910 ALLOCATE( & 911 & zsurfm(imaxifp,imaxjfp,isurf) & 912 & ) 913 914 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 915 & surfdataqc%vdmean(:,:), zsurfm ) 916 917 ENDIF 558 918 559 919 ! Loop over observations 560 561 DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 562 563 iobs = jobs - sladatqc%nsurfup 564 565 IF ( kt /= sladatqc%mstp(jobs) ) THEN 566 920 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 921 922 iobs = jobs - surfdataqc%nsurfup 923 924 IF ( kt /= surfdataqc%mstp(jobs) ) THEN 925 567 926 IF(lwp) THEN 568 927 WRITE(numout,*) … … 574 933 WRITE(numout,*) ' Record = ', jobs, & 575 934 & ' kt = ', kt, & 576 & ' mstp = ', s ladatqc%mstp(jobs), &577 & ' ntyp = ', s ladatqc%ntyp(jobs)935 & ' mstp = ', surfdataqc%mstp(jobs), & 936 & ' ntyp = ', surfdataqc%ntyp(jobs) 578 937 ENDIF 579 CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 580 938 CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 939 940 ENDIF 941 942 zlam = surfdataqc%rlam(jobs) 943 zphi = surfdataqc%rphi(jobs) 944 945 IF ( ldnightav .AND. idayend == 0 ) THEN 946 ! Night-time averaged data 947 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 948 ELSE 949 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 950 ENDIF 951 952 IF ( k2dint <= 4 ) THEN 953 954 ! Get weights to interpolate the model value to the observation point 955 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 956 & zglam(:,:,iobs), zgphi(:,:,iobs), & 957 & zmask(:,:,iobs), zweig, zobsmask ) 958 959 ! Interpolate the model value to the observation point 960 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 961 962 ELSE 963 964 ! Get weights to average the model SLA to the observation footprint 965 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 966 & zglam(:,:,iobs), zgphi(:,:,iobs), & 967 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 968 & zmask(:,:,iobs), plamscl, pphiscl, & 969 & lindegrees, zweig, zobsmask ) 970 971 ! Average the model SST to the observation footprint 972 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 973 & zweig, zsurftmp(:,:,iobs), zext ) 974 975 ENDIF 976 977 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 978 ! ... Remove the MDT from the SSH at the observation point to get the SLA 979 surfdataqc%rext(jobs,1) = zext(1) 980 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 981 ELSE 982 surfdataqc%rmod(jobs,1) = zext(1) 581 983 ENDIF 582 984 583 zlam = sladatqc%rlam(jobs) 584 zphi = sladatqc%rphi(jobs) 585 586 ! Get weights to interpolate the model SSH to the observation point 587 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 588 & zglam(:,:,iobs), zgphi(:,:,iobs), & 589 & zmask(:,:,iobs), zweig, zobsmask ) 590 591 592 ! Interpolate the model SSH to the observation point 593 CALL obs_int_h2d( 1, 1, & 594 & zweig, zsshl(:,:,iobs), zext ) 595 596 sladatqc%rext(jobs,1) = zext(1) 597 ! ... Remove the MDT at the observation point 598 sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 985 IF ( zext(1) == obfillflt ) THEN 986 ! If the observation value is a fill value, set QC flag to bad 987 surfdataqc%nqc(jobs) = 4 988 ENDIF 599 989 600 990 END DO … … 602 992 ! Deallocate the data for interpolation 603 993 DEALLOCATE( & 994 & zweig, & 604 995 & igrdi, & 605 996 & igrdj, & … … 607 998 & zgphi, & 608 999 & zmask, & 609 & zsshl & 1000 & zsurf, & 1001 & zsurftmp, & 1002 & zglamf, & 1003 & zgphif, & 1004 & igrdip1,& 1005 & igrdjp1 & 610 1006 & ) 611 1007 612 sladatqc%nsurfup = sladatqc%nsurfup + isla 613 614 END SUBROUTINE obs_sla_opt 615 616 SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 617 & psstn, psstmask, k2dint, ld_nightav ) 618 !!----------------------------------------------------------------------- 619 !! 620 !! *** ROUTINE obs_sst_opt *** 621 !! 622 !! ** Purpose : Compute the model counterpart of surface temperature 623 !! data by interpolating from the model grid to the 624 !! observation point. 625 !! 626 !! ** Method : Linearly interpolate to each observation point using 627 !! the model values at the corners of the surrounding grid box. 628 !! 629 !! The now model SST is first computed at the obs (lon, lat) point. 630 !! 631 !! Several horizontal interpolation schemes are available: 632 !! - distance-weighted (great circle) (k2dint = 0) 633 !! - distance-weighted (small angle) (k2dint = 1) 634 !! - bilinear (geographical grid) (k2dint = 2) 635 !! - bilinear (quadrilateral grid) (k2dint = 3) 636 !! - polynomial (quadrilateral grid) (k2dint = 4) 637 !! 638 !! 639 !! ** Action : 640 !! 641 !! History : 642 !! ! 07-07 (S. Ricci ) : Original 643 !! 644 !!----------------------------------------------------------------------- 645 646 !! * Modules used 647 USE obs_surf_def ! Definition of storage space for surface observations 648 USE sbcdcy 649 650 IMPLICIT NONE 651 652 !! * Arguments 653 TYPE(obs_surf), INTENT(INOUT) :: & 654 & sstdatqc ! Subset of surface data not failing screening 655 INTEGER, INTENT(IN) :: kt ! Time step 656 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 657 INTEGER, INTENT(IN) :: kpj 658 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 659 ! (kit000-1 = restart time) 660 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 661 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 662 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 663 & psstn, & ! Model SST field 664 & psstmask ! Land-sea mask 665 666 !! * Local declarations 667 INTEGER :: ji 668 INTEGER :: jj 669 INTEGER :: jobs 670 INTEGER :: inrc 671 INTEGER :: isst 672 INTEGER :: iobs 673 INTEGER :: idayend 674 REAL(KIND=wp) :: zlam 675 REAL(KIND=wp) :: zphi 676 REAL(KIND=wp) :: zext(1), zobsmask(1) 677 REAL(KIND=wp) :: zdaystp 678 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 679 & icount_sstnight, & 680 & imask_night 681 REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 682 & zintmp, & 683 & zouttmp, & 684 & zmeanday ! to compute model sst in region of 24h daylight (pole) 685 REAL(kind=wp), DIMENSION(2,2,1) :: & 686 & zweig 687 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 688 & zmask, & 689 & zsstl, & 690 & zsstm, & 691 & zglam, & 692 & zgphi 693 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 694 & igrdi, & 695 & igrdj 696 LOGICAL, INTENT(IN) :: ld_nightav 697 698 !----------------------------------------------------------------------- 699 ! Local initialization 700 !----------------------------------------------------------------------- 701 ! ... Record and data counters 702 inrc = kt - kit000 + 2 703 isst = sstdatqc%nsstp(inrc) 704 705 IF ( ld_nightav ) THEN 706 707 ! Initialize array for night mean 708 709 IF ( kt .EQ. 0 ) THEN 710 ALLOCATE ( icount_sstnight(kpi,kpj) ) 711 ALLOCATE ( imask_night(kpi,kpj) ) 712 ALLOCATE ( zintmp(kpi,kpj) ) 713 ALLOCATE ( zouttmp(kpi,kpj) ) 714 ALLOCATE ( zmeanday(kpi,kpj) ) 715 nday_qsr = -1 ! initialisation flag for nbc_dcy 716 ENDIF 717 718 ! Initialize daily mean for first timestep 719 idayend = MOD( kt - kit000 + 1, kdaystp ) 720 721 ! Added kt == 0 test to catch restart case 722 IF ( idayend == 1 .OR. kt == 0) THEN 723 IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 724 DO jj = 1, jpj 725 DO ji = 1, jpi 726 sstdatqc%vdmean(ji,jj) = 0.0 727 zmeanday(ji,jj) = 0.0 728 icount_sstnight(ji,jj) = 0 729 END DO 730 END DO 731 ENDIF 732 733 zintmp(:,:) = 0.0 734 zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 735 imask_night(:,:) = INT( zouttmp(:,:) ) 736 737 DO jj = 1, jpj 738 DO ji = 1, jpi 739 ! Increment the temperature field for computing night mean and counter 740 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 741 & + psstn(ji,jj)*imask_night(ji,jj) 742 zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj) 743 icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 744 END DO 745 END DO 746 747 ! Compute the daily mean at the end of day 748 749 zdaystp = 1.0 / REAL( kdaystp ) 750 751 IF ( idayend == 0 ) THEN 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 ! Test if "no night" point 755 IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 756 sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 757 & / icount_sstnight(ji,jj) 758 ELSE 759 sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 760 ENDIF 761 END DO 762 END DO 763 ENDIF 764 765 ENDIF 766 767 ! Get the data for interpolation 768 769 ALLOCATE( & 770 & igrdi(2,2,isst), & 771 & igrdj(2,2,isst), & 772 & zglam(2,2,isst), & 773 & zgphi(2,2,isst), & 774 & zmask(2,2,isst), & 775 & zsstl(2,2,isst) & 776 & ) 777 778 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 779 iobs = jobs - sstdatqc%nsurfup 780 igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 781 igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 782 igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 783 igrdj(1,2,iobs) = sstdatqc%mj(jobs) 784 igrdi(2,1,iobs) = sstdatqc%mi(jobs) 785 igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 786 igrdi(2,2,iobs) = sstdatqc%mi(jobs) 787 igrdj(2,2,iobs) = sstdatqc%mj(jobs) 788 END DO 789 790 CALL obs_int_comm_2d( 2, 2, isst, & 791 & igrdi, igrdj, glamt, zglam ) 792 CALL obs_int_comm_2d( 2, 2, isst, & 793 & igrdi, igrdj, gphit, zgphi ) 794 CALL obs_int_comm_2d( 2, 2, isst, & 795 & igrdi, igrdj, psstmask, zmask ) 796 CALL obs_int_comm_2d( 2, 2, isst, & 797 & igrdi, igrdj, psstn, zsstl ) 798 799 ! At the end of the day get interpolated means 800 IF ( idayend == 0 .AND. ld_nightav ) THEN 801 802 ALLOCATE( & 803 & zsstm(2,2,isst) & 804 & ) 805 806 CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 807 & sstdatqc%vdmean(:,:), zsstm ) 808 809 ENDIF 810 811 ! Loop over observations 812 813 DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 814 815 iobs = jobs - sstdatqc%nsurfup 816 817 IF ( kt /= sstdatqc%mstp(jobs) ) THEN 818 819 IF(lwp) THEN 820 WRITE(numout,*) 821 WRITE(numout,*) ' E R R O R : Observation', & 822 & ' time step is not consistent with the', & 823 & ' model time step' 824 WRITE(numout,*) ' =========' 825 WRITE(numout,*) 826 WRITE(numout,*) ' Record = ', jobs, & 827 & ' kt = ', kt, & 828 & ' mstp = ', sstdatqc%mstp(jobs), & 829 & ' ntyp = ', sstdatqc%ntyp(jobs) 830 ENDIF 831 CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 832 833 ENDIF 834 835 zlam = sstdatqc%rlam(jobs) 836 zphi = sstdatqc%rphi(jobs) 837 838 ! Get weights to interpolate the model SST to the observation point 839 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 840 & zglam(:,:,iobs), zgphi(:,:,iobs), & 841 & zmask(:,:,iobs), zweig, zobsmask ) 842 843 ! Interpolate the model SST to the observation point 844 845 IF ( ld_nightav ) THEN 846 847 IF ( idayend == 0 ) THEN 848 ! Daily averaged/diurnal cycle of SST data 849 CALL obs_int_h2d( 1, 1, & 850 & zweig, zsstm(:,:,iobs), zext ) 851 ELSE 852 CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // & 853 & ' number of night SST data should' // & 854 & ' only occur at the end of a given day' ) 855 ENDIF 856 857 ELSE 858 859 CALL obs_int_h2d( 1, 1, & 860 & zweig, zsstl(:,:,iobs), zext ) 861 862 ENDIF 863 sstdatqc%rmod(jobs,1) = zext(1) 864 865 END DO 866 867 ! Deallocate the data for interpolation 868 DEALLOCATE( & 869 & igrdi, & 870 & igrdj, & 871 & zglam, & 872 & zgphi, & 873 & zmask, & 874 & zsstl & 875 & ) 876 877 ! At the end of the day also get interpolated means 878 IF ( idayend == 0 .AND. ld_nightav ) THEN 1008 ! At the end of the day also deallocate night-time mean array 1009 IF ( idayend == 0 .AND. ldnightav ) THEN 879 1010 DEALLOCATE( & 880 & zs stm &1011 & zsurfm & 881 1012 & ) 882 1013 ENDIF 883 884 sstdatqc%nsurfup = sstdatqc%nsurfup + isst 885 886 END SUBROUTINE obs_sst_opt 887 888 SUBROUTINE obs_sss_opt 889 !!----------------------------------------------------------------------- 890 !! 891 !! *** ROUTINE obs_sss_opt *** 892 !! 893 !! ** Purpose : Compute the model counterpart of sea surface salinity 894 !! data by interpolating from the model grid to the 895 !! observation point. 896 !! 897 !! ** Method : 898 !! 899 !! ** Action : 900 !! 901 !! History : 902 !! ! ??-?? 903 !!----------------------------------------------------------------------- 904 905 IMPLICIT NONE 906 907 END SUBROUTINE obs_sss_opt 908 909 SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, & 910 & pseaicen, pseaicemask, k2dint ) 911 912 !!----------------------------------------------------------------------- 913 !! 914 !! *** ROUTINE obs_seaice_opt *** 915 !! 916 !! ** Purpose : Compute the model counterpart of surface temperature 917 !! data by interpolating from the model grid to the 918 !! observation point. 919 !! 920 !! ** Method : Linearly interpolate to each observation point using 921 !! the model values at the corners of the surrounding grid box. 922 !! 923 !! The now model sea ice is first computed at the obs (lon, lat) point. 924 !! 925 !! Several horizontal interpolation schemes are available: 926 !! - distance-weighted (great circle) (k2dint = 0) 927 !! - distance-weighted (small angle) (k2dint = 1) 928 !! - bilinear (geographical grid) (k2dint = 2) 929 !! - bilinear (quadrilateral grid) (k2dint = 3) 930 !! - polynomial (quadrilateral grid) (k2dint = 4) 931 !! 932 !! 933 !! ** Action : 934 !! 935 !! History : 936 !! ! 07-07 (S. Ricci ) : Original 937 !! 938 !!----------------------------------------------------------------------- 939 940 !! * Modules used 941 USE obs_surf_def ! Definition of storage space for surface observations 942 943 IMPLICIT NONE 944 945 !! * Arguments 946 TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening 947 INTEGER, INTENT(IN) :: kt ! Time step 948 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 949 INTEGER, INTENT(IN) :: kpj 950 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 951 ! (kit000-1 = restart time) 952 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 953 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 954 & pseaicen, & ! Model sea ice field 955 & pseaicemask ! Land-sea mask 956 957 !! * Local declarations 958 INTEGER :: ji 959 INTEGER :: jj 960 INTEGER :: jobs 961 INTEGER :: inrc 962 INTEGER :: iseaice 963 INTEGER :: iobs 964 965 REAL(KIND=wp) :: zlam 966 REAL(KIND=wp) :: zphi 967 REAL(KIND=wp) :: zext(1), zobsmask(1) 968 REAL(kind=wp), DIMENSION(2,2,1) :: & 969 & zweig 970 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 971 & zmask, & 972 & zseaicel, & 973 & zglam, & 974 & zgphi 975 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 976 & igrdi, & 977 & igrdj 978 979 !------------------------------------------------------------------------ 980 ! Local initialization 981 !------------------------------------------------------------------------ 982 ! ... Record and data counters 983 inrc = kt - kit000 + 2 984 iseaice = seaicedatqc%nsstp(inrc) 985 986 ! Get the data for interpolation 987 988 ALLOCATE( & 989 & igrdi(2,2,iseaice), & 990 & igrdj(2,2,iseaice), & 991 & zglam(2,2,iseaice), & 992 & zgphi(2,2,iseaice), & 993 & zmask(2,2,iseaice), & 994 & zseaicel(2,2,iseaice) & 995 & ) 996 997 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 998 iobs = jobs - seaicedatqc%nsurfup 999 igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1 1000 igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1 1001 igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1 1002 igrdj(1,2,iobs) = seaicedatqc%mj(jobs) 1003 igrdi(2,1,iobs) = seaicedatqc%mi(jobs) 1004 igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1 1005 igrdi(2,2,iobs) = seaicedatqc%mi(jobs) 1006 igrdj(2,2,iobs) = seaicedatqc%mj(jobs) 1007 END DO 1008 1009 CALL obs_int_comm_2d( 2, 2, iseaice, & 1010 & igrdi, igrdj, glamt, zglam ) 1011 CALL obs_int_comm_2d( 2, 2, iseaice, & 1012 & igrdi, igrdj, gphit, zgphi ) 1013 CALL obs_int_comm_2d( 2, 2, iseaice, & 1014 & igrdi, igrdj, pseaicemask, zmask ) 1015 CALL obs_int_comm_2d( 2, 2, iseaice, & 1016 & igrdi, igrdj, pseaicen, zseaicel ) 1017 1018 DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 1019 1020 iobs = jobs - seaicedatqc%nsurfup 1021 1022 IF ( kt /= seaicedatqc%mstp(jobs) ) THEN 1023 1024 IF(lwp) THEN 1025 WRITE(numout,*) 1026 WRITE(numout,*) ' E R R O R : Observation', & 1027 & ' time step is not consistent with the', & 1028 & ' model time step' 1029 WRITE(numout,*) ' =========' 1030 WRITE(numout,*) 1031 WRITE(numout,*) ' Record = ', jobs, & 1032 & ' kt = ', kt, & 1033 & ' mstp = ', seaicedatqc%mstp(jobs), & 1034 & ' ntyp = ', seaicedatqc%ntyp(jobs) 1035 ENDIF 1036 CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' ) 1037 1038 ENDIF 1039 1040 zlam = seaicedatqc%rlam(jobs) 1041 zphi = seaicedatqc%rphi(jobs) 1042 1043 ! Get weights to interpolate the model sea ice to the observation point 1044 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 1045 & zglam(:,:,iobs), zgphi(:,:,iobs), & 1046 & zmask(:,:,iobs), zweig, zobsmask ) 1047 1048 ! ... Interpolate the model sea ice to the observation point 1049 CALL obs_int_h2d( 1, 1, & 1050 & zweig, zseaicel(:,:,iobs), zext ) 1051 1052 seaicedatqc%rmod(jobs,1) = zext(1) 1053 1054 END DO 1055 1056 ! Deallocate the data for interpolation 1057 DEALLOCATE( & 1058 & igrdi, & 1059 & igrdj, & 1060 & zglam, & 1061 & zgphi, & 1062 & zmask, & 1063 & zseaicel & 1064 & ) 1065 1066 seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice 1067 1068 END SUBROUTINE obs_seaice_opt 1069 1070 SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 1071 & pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, & 1072 & ld_dailyav ) 1073 !!----------------------------------------------------------------------- 1074 !! 1075 !! *** ROUTINE obs_vel_opt *** 1076 !! 1077 !! ** Purpose : Compute the model counterpart of velocity profile 1078 !! data by interpolating from the model grid to the 1079 !! observation point. 1080 !! 1081 !! ** Method : Linearly interpolate zonal and meridional components of velocity 1082 !! to each observation point using the model values at the corners of 1083 !! the surrounding grid box. The model velocity components are on a 1084 !! staggered C- grid. 1085 !! 1086 !! For velocity data from the TAO array, the model equivalent is 1087 !! a daily mean velocity field. So, we first compute 1088 !! the mean, then interpolate only at the end of the day. 1089 !! 1090 !! ** Action : 1091 !! 1092 !! History : 1093 !! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles 1094 !! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles 1095 !!----------------------------------------------------------------------- 1096 1097 !! * Modules used 1098 USE obs_profiles_def ! Definition of storage space for profile obs. 1099 1100 IMPLICIT NONE 1101 1102 !! * Arguments 1103 TYPE(obs_prof), INTENT(INOUT) :: & 1104 & prodatqc ! Subset of profile data not failing screening 1105 INTEGER, INTENT(IN) :: kt ! Time step 1106 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 1107 INTEGER, INTENT(IN) :: kpj 1108 INTEGER, INTENT(IN) :: kpk 1109 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 1110 ! (kit000-1 = restart time) 1111 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 1112 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 1113 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 1114 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 1115 & pun, & ! Model zonal component of velocity 1116 & pvn, & ! Model meridional component of velocity 1117 & pumask, & ! Land-sea mask 1118 & pvmask ! Land-sea mask 1119 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 1120 & pgdept ! Model array of depth levels 1121 LOGICAL, INTENT(IN) :: ld_dailyav 1122 1123 !! * Local declarations 1124 INTEGER :: ji 1125 INTEGER :: jj 1126 INTEGER :: jk 1127 INTEGER :: jobs 1128 INTEGER :: inrc 1129 INTEGER :: ipro 1130 INTEGER :: idayend 1131 INTEGER :: ista 1132 INTEGER :: iend 1133 INTEGER :: iobs 1134 INTEGER, DIMENSION(imaxavtypes) :: & 1135 & idailyavtypes 1136 REAL(KIND=wp) :: zlam 1137 REAL(KIND=wp) :: zphi 1138 REAL(KIND=wp) :: zdaystp 1139 REAL(KIND=wp), DIMENSION(kpk) :: & 1140 & zobsmasku, & 1141 & zobsmaskv, & 1142 & zobsmask, & 1143 & zobsk, & 1144 & zobs2k 1145 REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 1146 & zweigu,zweigv 1147 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 1148 & zumask, zvmask, & 1149 & zintu, & 1150 & zintv, & 1151 & zinmu, & 1152 & zinmv 1153 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 1154 & zglamu, zglamv, & 1155 & zgphiu, zgphiv 1156 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1157 & igrdiu, & 1158 & igrdju, & 1159 & igrdiv, & 1160 & igrdjv 1161 1162 !------------------------------------------------------------------------ 1163 ! Local initialization 1164 !------------------------------------------------------------------------ 1165 ! ... Record and data counters 1166 inrc = kt - kit000 + 2 1167 ipro = prodatqc%npstp(inrc) 1168 1169 ! Initialize daily mean for first timestep 1170 idayend = MOD( kt - kit000 + 1, kdaystp ) 1171 1172 ! Added kt == 0 test to catch restart case 1173 IF ( idayend == 1 .OR. kt == 0) THEN 1174 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 1175 prodatqc%vdmean(:,:,:,1) = 0.0 1176 prodatqc%vdmean(:,:,:,2) = 0.0 1177 ENDIF 1178 1179 ! Increment the zonal velocity field for computing daily mean 1180 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:) 1181 ! Increment the meridional velocity field for computing daily mean 1182 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:) 1183 1184 ! Compute the daily mean at the end of day 1185 zdaystp = 1.0 / REAL( kdaystp ) 1186 IF ( idayend == 0 ) THEN 1187 prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp 1188 prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp 1189 ENDIF 1190 1191 ! Get the data for interpolation 1192 ALLOCATE( & 1193 & igrdiu(2,2,ipro), & 1194 & igrdju(2,2,ipro), & 1195 & igrdiv(2,2,ipro), & 1196 & igrdjv(2,2,ipro), & 1197 & zglamu(2,2,ipro), zglamv(2,2,ipro), & 1198 & zgphiu(2,2,ipro), zgphiv(2,2,ipro), & 1199 & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), & 1200 & zintu(2,2,kpk,ipro), & 1201 & zintv(2,2,kpk,ipro) & 1202 & ) 1203 1204 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1205 iobs = jobs - prodatqc%nprofup 1206 igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1 1207 igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1 1208 igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1 1209 igrdju(1,2,iobs) = prodatqc%mj(jobs,1) 1210 igrdiu(2,1,iobs) = prodatqc%mi(jobs,1) 1211 igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1 1212 igrdiu(2,2,iobs) = prodatqc%mi(jobs,1) 1213 igrdju(2,2,iobs) = prodatqc%mj(jobs,1) 1214 igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1 1215 igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1 1216 igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1 1217 igrdjv(1,2,iobs) = prodatqc%mj(jobs,2) 1218 igrdiv(2,1,iobs) = prodatqc%mi(jobs,2) 1219 igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1 1220 igrdiv(2,2,iobs) = prodatqc%mi(jobs,2) 1221 igrdjv(2,2,iobs) = prodatqc%mj(jobs,2) 1222 END DO 1223 1224 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu ) 1225 CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu ) 1226 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask ) 1227 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu ) 1228 1229 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv ) 1230 CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv ) 1231 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask ) 1232 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv ) 1233 1234 ! At the end of the day also get interpolated means 1235 IF ( idayend == 0 ) THEN 1236 1237 ALLOCATE( & 1238 & zinmu(2,2,kpk,ipro), & 1239 & zinmv(2,2,kpk,ipro) & 1240 & ) 1241 1242 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, & 1243 & prodatqc%vdmean(:,:,:,1), zinmu ) 1244 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, & 1245 & prodatqc%vdmean(:,:,:,2), zinmv ) 1246 1247 ENDIF 1248 1249 ! loop over observations 1250 1251 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 1252 1253 iobs = jobs - prodatqc%nprofup 1254 1255 IF ( kt /= prodatqc%mstp(jobs) ) THEN 1256 1257 IF(lwp) THEN 1258 WRITE(numout,*) 1259 WRITE(numout,*) ' E R R O R : Observation', & 1260 & ' time step is not consistent with the', & 1261 & ' model time step' 1262 WRITE(numout,*) ' =========' 1263 WRITE(numout,*) 1264 WRITE(numout,*) ' Record = ', jobs, & 1265 & ' kt = ', kt, & 1266 & ' mstp = ', prodatqc%mstp(jobs), & 1267 & ' ntyp = ', prodatqc%ntyp(jobs) 1268 ENDIF 1269 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 1270 ENDIF 1271 1272 zlam = prodatqc%rlam(jobs) 1273 zphi = prodatqc%rphi(jobs) 1274 1275 ! Initialize observation masks 1276 1277 zobsmasku(:) = 0.0 1278 zobsmaskv(:) = 0.0 1279 1280 ! Horizontal weights and vertical mask 1281 1282 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1283 1284 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1285 & zglamu(:,:,iobs), zgphiu(:,:,iobs), & 1286 & zumask(:,:,:,iobs), zweigu, zobsmasku ) 1287 1288 ENDIF 1289 1290 1291 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1292 1293 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1294 & zglamv(:,:,iobs), zgphiv(:,:,iobs), & 1295 & zvmask(:,:,:,iobs), zweigv, zobsmasku ) 1296 1297 ENDIF 1298 1299 ! Ensure that the vertical mask on u and v are consistent. 1300 1301 zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) ) 1302 1303 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 1304 1305 zobsk(:) = obfillflt 1306 1307 IF ( ld_dailyav ) THEN 1308 1309 IF ( idayend == 0 ) THEN 1310 1311 ! Daily averaged data 1312 1313 CALL obs_int_h2d( kpk, kpk, & 1314 & zweigu, zinmu(:,:,:,iobs), zobsk ) 1315 1316 1317 ELSE 1318 1319 CALL ctl_stop( ' A nonzero' // & 1320 & ' number of U profile data should' // & 1321 & ' only occur at the end of a given day' ) 1322 1323 ENDIF 1324 1325 ELSE 1326 1327 ! Point data 1328 1329 CALL obs_int_h2d( kpk, kpk, & 1330 & zweigu, zintu(:,:,:,iobs), zobsk ) 1331 1332 ENDIF 1333 1334 !------------------------------------------------------------- 1335 ! Compute vertical second-derivative of the interpolating 1336 ! polynomial at obs points 1337 !------------------------------------------------------------- 1338 1339 IF ( k1dint == 1 ) THEN 1340 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 1341 & pgdept, zobsmask ) 1342 ENDIF 1343 1344 !----------------------------------------------------------------- 1345 ! Vertical interpolation to the observation point 1346 !----------------------------------------------------------------- 1347 ista = prodatqc%npvsta(jobs,1) 1348 iend = prodatqc%npvend(jobs,1) 1349 CALL obs_int_z1d( kpk, & 1350 & prodatqc%var(1)%mvk(ista:iend), & 1351 & k1dint, iend - ista + 1, & 1352 & prodatqc%var(1)%vdep(ista:iend), & 1353 & zobsk, zobs2k, & 1354 & prodatqc%var(1)%vmod(ista:iend), & 1355 & pgdept, zobsmask ) 1356 1357 ENDIF 1358 1359 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 1360 1361 zobsk(:) = obfillflt 1362 1363 IF ( ld_dailyav ) THEN 1364 1365 IF ( idayend == 0 ) THEN 1366 1367 ! Daily averaged data 1368 1369 CALL obs_int_h2d( kpk, kpk, & 1370 & zweigv, zinmv(:,:,:,iobs), zobsk ) 1371 1372 ELSE 1373 1374 CALL ctl_stop( ' A nonzero' // & 1375 & ' number of V profile data should' // & 1376 & ' only occur at the end of a given day' ) 1377 1378 ENDIF 1379 1380 ELSE 1381 1382 ! Point data 1383 1384 CALL obs_int_h2d( kpk, kpk, & 1385 & zweigv, zintv(:,:,:,iobs), zobsk ) 1386 1387 ENDIF 1388 1389 1390 !------------------------------------------------------------- 1391 ! Compute vertical second-derivative of the interpolating 1392 ! polynomial at obs points 1393 !------------------------------------------------------------- 1394 1395 IF ( k1dint == 1 ) THEN 1396 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 1397 & pgdept, zobsmask ) 1398 ENDIF 1399 1400 !---------------------------------------------------------------- 1401 ! Vertical interpolation to the observation point 1402 !---------------------------------------------------------------- 1403 ista = prodatqc%npvsta(jobs,2) 1404 iend = prodatqc%npvend(jobs,2) 1405 CALL obs_int_z1d( kpk, & 1406 & prodatqc%var(2)%mvk(ista:iend),& 1407 & k1dint, iend - ista + 1, & 1408 & prodatqc%var(2)%vdep(ista:iend),& 1409 & zobsk, zobs2k, & 1410 & prodatqc%var(2)%vmod(ista:iend),& 1411 & pgdept, zobsmask ) 1412 1413 ENDIF 1414 1415 END DO 1416 1417 ! Deallocate the data for interpolation 1418 DEALLOCATE( & 1419 & igrdiu, & 1420 & igrdju, & 1421 & igrdiv, & 1422 & igrdjv, & 1423 & zglamu, zglamv, & 1424 & zgphiu, zgphiv, & 1425 & zumask, zvmask, & 1426 & zintu, & 1427 & zintv & 1428 & ) 1429 ! At the end of the day also get interpolated means 1430 IF ( idayend == 0 ) THEN 1431 DEALLOCATE( & 1432 & zinmu, & 1433 & zinmv & 1434 & ) 1435 ENDIF 1436 1437 prodatqc%nprofup = prodatqc%nprofup + ipro 1438 1439 END SUBROUTINE obs_vel_opt 1014 1015 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1016 1017 END SUBROUTINE obs_surf_opt 1440 1018 1441 1019 END MODULE obs_oper 1442
Note: See TracChangeset
for help on using the changeset viewer.
