Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r4990 r6140 25 25 USE netcdf ! NetCDF library 26 26 USE obs_oper ! Observation operators 27 USE obs_prof_io ! Profile files I/O (non-FB files)28 27 USE lib_mpp ! For ctl_warn/stop 28 USE obs_fbm ! Feedback routines 29 29 30 30 IMPLICIT NONE … … 33 33 PRIVATE 34 34 35 PUBLIC obs_rea_pro _dri! Read the profile observations35 PUBLIC obs_rea_prof ! Read the profile observations 36 36 37 37 !!---------------------------------------------------------------------- … … 42 42 43 43 CONTAINS 44 45 SUBROUTINE obs_rea_pro_dri( kformat, & 46 & profdata, knumfiles, cfilenames, & 47 & kvars, kextr, kstp, ddobsini, ddobsend, & 48 & ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, & 49 & ldmod, kdailyavtypes ) 44 45 SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, & 46 & kvars, kextr, kstp, ddobsini, ddobsend, & 47 & ldvar1, ldvar2, ldignmis, ldsatt, & 48 & ldmod, kdailyavtypes ) 50 49 !!--------------------------------------------------------------------- 51 50 !! 52 !! *** ROUTINE obs_rea_pro _dri***51 !! *** ROUTINE obs_rea_prof *** 53 52 !! 54 53 !! ** Purpose : Read from file the profile observations 55 54 !! 56 !! ** Method : Depending on kformat either ENACT, CORIOLIS or57 !! feedback data files are read55 !! ** Method : Read feedback data in and transform to NEMO internal 56 !! profile data structure 58 57 !! 59 58 !! ** Action : … … 63 62 !! History : 64 63 !! ! : 2009-09 (K. Mogensen) : New merged version of old routines 64 !! ! : 2015-08 (M. Martin) : Merged profile and velocity routines 65 65 !!---------------------------------------------------------------------- 66 !! * Modules used 67 66 68 67 !! * Arguments 69 INTEGER :: kformat ! Format of input data 70 ! ! 1: ENACT 71 ! ! 2: Coriolis 72 TYPE(obs_prof), INTENT(OUT) :: profdata ! Profile data to be read 73 INTEGER, INTENT(IN) :: knumfiles ! Number of files to read in 68 TYPE(obs_prof), INTENT(OUT) :: & 69 & profdata ! Profile data to be read 70 INTEGER, INTENT(IN) :: knumfiles ! Number of files to read 74 71 CHARACTER(LEN=128), INTENT(IN) :: & 75 & c filenames(knumfiles)! File names to read in72 & cdfilenames(knumfiles) ! File names to read in 76 73 INTEGER, INTENT(IN) :: kvars ! Number of variables in profdata 77 INTEGER, INTENT(IN) :: kextr ! Number of extra fields for each var in profdata 78 INTEGER, INTENT(IN) :: kstp ! Ocean time-step index 79 LOGICAL, INTENT(IN) :: ldt3d ! Observed variables switches 80 LOGICAL, INTENT(IN) :: lds3d 81 LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files 82 LOGICAL, INTENT(IN) :: ldsatt ! Compute salinity at all temperature points 83 LOGICAL, INTENT(IN) :: ldavtimset ! Correct time for daily averaged data 84 LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data 85 REAL(KIND=dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS 86 REAL(KIND=dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS 74 INTEGER, INTENT(IN) :: kextr ! Number of extra fields for each var 75 INTEGER, INTENT(IN) :: kstp ! Ocean time-step index 76 LOGICAL, INTENT(IN) :: ldvar1 ! Observed variables switches 77 LOGICAL, INTENT(IN) :: ldvar2 78 LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files 79 LOGICAL, INTENT(IN) :: ldsatt ! Compute salinity at all temperature points 80 LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data 81 REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS 82 REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS 87 83 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 88 & kdailyavtypes 84 & kdailyavtypes ! Types of daily average observations 89 85 90 86 !! * Local declarations 91 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_pro_dri' 87 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 88 CHARACTER(len=8) :: clrefdate 89 CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 92 90 INTEGER :: jvar 93 91 INTEGER :: ji … … 105 103 INTEGER :: imin 106 104 INTEGER :: isec 105 INTEGER :: iprof 106 INTEGER :: iproftot 107 INTEGER :: ivar1t0 108 INTEGER :: ivar2t0 109 INTEGER :: ivar1t 110 INTEGER :: ivar2t 111 INTEGER :: ip3dt 112 INTEGER :: ios 113 INTEGER :: ioserrcount 114 INTEGER :: ivar1tmpp 115 INTEGER :: ivar2tmpp 116 INTEGER :: ip3dtmpp 117 INTEGER :: itype 107 118 INTEGER, DIMENSION(knumfiles) :: & 108 119 & irefdate 109 120 INTEGER, DIMENSION(ntyp1770+1) :: & 110 & itypt, & 111 & ityptmpp, & 112 & ityps, & 113 & itypsmpp 114 INTEGER :: it3dtmpp 115 INTEGER :: is3dtmpp 116 INTEGER :: ip3dtmpp 121 & itypvar1, & 122 & itypvar1mpp, & 123 & itypvar2, & 124 & itypvar2mpp 117 125 INTEGER, DIMENSION(:), ALLOCATABLE :: & 118 & iobsi, & 119 & iobsj, & 120 & iproc, & 126 & iobsi1, & 127 & iobsj1, & 128 & iproc1, & 129 & iobsi2, & 130 & iobsj2, & 131 & iproc2, & 121 132 & iindx, & 122 133 & ifileidx, & 123 134 & iprofidx 124 INTEGER :: itype125 135 INTEGER, DIMENSION(imaxavtypes) :: & 126 136 & idailyavtypes 137 INTEGER, DIMENSION(kvars) :: & 138 & iv3dt 127 139 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 128 140 & zphi, & 129 141 & zlam 130 real(wp), DIMENSION(:), ALLOCATABLE :: &142 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 131 143 & zdat 144 REAL(wp), DIMENSION(knumfiles) :: & 145 & djulini, & 146 & djulend 132 147 LOGICAL :: llvalprof 148 LOGICAL :: lldavtimset 133 149 TYPE(obfbdata), POINTER, DIMENSION(:) :: & 134 150 & inpfiles 135 real(wp), DIMENSION(knumfiles) :: & 136 & djulini, & 137 & djulend 138 INTEGER :: iprof 139 INTEGER :: iproftot 140 INTEGER :: it3dt0 141 INTEGER :: is3dt0 142 INTEGER :: it3dt 143 INTEGER :: is3dt 144 INTEGER :: ip3dt 145 INTEGER :: ios 146 INTEGER :: ioserrcount 147 INTEGER, DIMENSION(kvars) :: & 148 & iv3dt 149 CHARACTER(len=8) :: cl_refdate 150 151 151 152 ! Local initialization 152 153 iprof = 0 153 i t3dt0 = 0154 i s3dt0 = 0154 ivar1t0 = 0 155 ivar2t0 = 0 155 156 ip3dt = 0 156 157 157 158 ! Daily average types 159 lldavtimset = .FALSE. 158 160 IF ( PRESENT(kdailyavtypes) ) THEN 159 161 idailyavtypes(:) = kdailyavtypes(:) 162 IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE. 160 163 ELSE 161 164 idailyavtypes(:) = -1 … … 163 166 164 167 !----------------------------------------------------------------------- 165 ! Check data the model part is just with feedback data files166 !-----------------------------------------------------------------------167 IF ( ldmod .AND. ( kformat /= 0 ) ) THEN168 CALL ctl_stop( 'Model can only be read from feedback data' )169 RETURN170 ENDIF171 172 !-----------------------------------------------------------------------173 168 ! Count the number of files needed and allocate the obfbdata type 174 169 !----------------------------------------------------------------------- 175 170 176 171 inobf = knumfiles 177 172 178 173 ALLOCATE( inpfiles(inobf) ) 179 174 180 175 prof_files : DO jj = 1, inobf 181 176 182 177 !--------------------------------------------------------------------- 183 178 ! Prints … … 186 181 WRITE(numout,*) 187 182 WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', & 188 & TRIM( TRIM( c filenames(jj) ) )183 & TRIM( TRIM( cdfilenames(jj) ) ) 189 184 WRITE(numout,*) ' ~~~~~~~~~~~~~~~' 190 185 WRITE(numout,*) … … 194 189 ! Initialization: Open file and get dimensions only 195 190 !--------------------------------------------------------------------- 196 197 iflag = nf90_open( TRIM( TRIM( cfilenames(jj)) ), nf90_nowrite, &191 192 iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, & 198 193 & i_file_id ) 199 194 200 195 IF ( iflag /= nf90_noerr ) THEN 201 196 202 197 IF ( ldignmis ) THEN 203 198 inpfiles(jj)%nobs = 0 204 CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj)) ) // &199 CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // & 205 200 & ' not found' ) 206 201 ELSE 207 CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj)) ) // &202 CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // & 208 203 & ' not found' ) 209 204 ENDIF 210 205 211 206 ELSE 212 207 213 208 !------------------------------------------------------------------ 214 ! Close the file since it is opened in read_ proffile209 ! Close the file since it is opened in read_obfbdata 215 210 !------------------------------------------------------------------ 216 211 217 212 iflag = nf90_close( i_file_id ) 218 213 … … 220 215 ! Read the profile file into inpfiles 221 216 !------------------------------------------------------------------ 222 IF ( kformat == 0 ) THEN 223 CALL init_obfbdata( inpfiles(jj) ) 224 IF(lwp) THEN 225 WRITE(numout,*) 226 WRITE(numout,*)'Reading from feedback file :', & 227 & TRIM( cfilenames(jj) ) 228 ENDIF 229 CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), & 230 & ldgrid = .TRUE. ) 231 IF ( inpfiles(jj)%nvar < 2 ) THEN 232 CALL ctl_stop( 'Feedback format error' ) 233 RETURN 234 ENDIF 235 IF ( TRIM(inpfiles(jj)%cname(1)) /= 'POTM' ) THEN 236 CALL ctl_stop( 'Feedback format error' ) 237 RETURN 238 ENDIF 239 IF ( TRIM(inpfiles(jj)%cname(2)) /= 'PSAL' ) THEN 240 CALL ctl_stop( 'Feedback format error' ) 241 RETURN 242 ENDIF 243 IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 244 CALL ctl_stop( 'Model not in input data' ) 245 RETURN 246 ENDIF 247 ELSEIF ( kformat == 1 ) THEN 248 CALL read_enactfile( TRIM( cfilenames(jj) ), inpfiles(jj), & 249 & numout, lwp, .TRUE. ) 250 ELSEIF ( kformat == 2 ) THEN 251 CALL read_coriofile( TRIM( cfilenames(jj) ), inpfiles(jj), & 252 & numout, lwp, .TRUE. ) 217 CALL init_obfbdata( inpfiles(jj) ) 218 CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), & 219 & ldgrid = .TRUE. ) 220 221 IF ( inpfiles(jj)%nvar < 2 ) THEN 222 CALL ctl_stop( 'Feedback format error: ', & 223 & ' less than 2 vars in profile file' ) 224 ENDIF 225 226 IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 227 CALL ctl_stop( 'Model not in input data' ) 228 ENDIF 229 230 IF ( jj == 1 ) THEN 231 ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 232 DO ji = 1, inpfiles(jj)%nvar 233 clvars(ji) = inpfiles(jj)%cname(ji) 234 END DO 253 235 ELSE 254 CALL ctl_stop( 'File format unknown' ) 255 ENDIF 256 236 DO ji = 1, inpfiles(jj)%nvar 237 IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN 238 CALL ctl_stop( 'Feedback file variables not consistent', & 239 & ' with previous files for this type' ) 240 ENDIF 241 END DO 242 ENDIF 243 257 244 !------------------------------------------------------------------ 258 245 ! Change longitude (-180,180) … … 272 259 ! Calculate the date (change eventually) 273 260 !------------------------------------------------------------------ 274 cl _refdate=inpfiles(jj)%cdjuldref(1:8)275 READ(cl _refdate,'(I8)') irefdate(jj)276 261 clrefdate=inpfiles(jj)%cdjuldref(1:8) 262 READ(clrefdate,'(I8)') irefdate(jj) 263 277 264 CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec ) 278 265 CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), & … … 283 270 284 271 ioserrcount=0 285 IF ( ldavtimset ) THEN 272 IF ( lldavtimset ) THEN 273 274 IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN 275 WRITE(numout,*)' Resetting time of daily averaged', & 276 & ' observations to the end of the day' 277 ENDIF 278 286 279 DO ji = 1, inpfiles(jj)%nobs 287 !288 ! for daily averaged data for example289 ! MRB data (itype==820) force the time290 ! to be the end of the day291 !292 280 READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype 293 281 900 IF ( ios /= 0 ) THEN 294 itype = 0 ! Set type to zero if there is a problem in the string conversion 295 ENDIF 296 IF ( ANY (idailyavtypes == itype ) ) THEN 297 inpfiles(jj)%ptim(ji) = & 298 & INT(inpfiles(jj)%ptim(ji)) + 1 299 ENDIF 282 ! Set type to zero if there is a problem in the string conversion 283 itype = 0 284 ENDIF 285 286 IF ( ANY ( idailyavtypes(:) == itype ) ) THEN 287 ! for daily averaged data force the time 288 ! to be the last time-step of the day, but still within the day. 289 IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN 290 inpfiles(jj)%ptim(ji) = & 291 & INT(inpfiles(jj)%ptim(ji)) + 0.9999 292 ELSE 293 inpfiles(jj)%ptim(ji) = & 294 & INT(inpfiles(jj)%ptim(ji)) - 0.0001 295 ENDIF 296 ENDIF 297 300 298 END DO 301 ENDIF 302 299 300 ENDIF 301 303 302 IF ( inpfiles(jj)%nobs > 0 ) THEN 304 inpfiles(jj)%iproc = -1305 inpfiles(jj)%iobsi = -1306 inpfiles(jj)%iobsj = -1303 inpfiles(jj)%iproc(:,:) = -1 304 inpfiles(jj)%iobsi(:,:) = -1 305 inpfiles(jj)%iobsj(:,:) = -1 307 306 ENDIF 308 307 inowin = 0 … … 318 317 ALLOCATE( zlam(inowin) ) 319 318 ALLOCATE( zphi(inowin) ) 320 ALLOCATE( iobsi(inowin) ) 321 ALLOCATE( iobsj(inowin) ) 322 ALLOCATE( iproc(inowin) ) 319 ALLOCATE( iobsi1(inowin) ) 320 ALLOCATE( iobsj1(inowin) ) 321 ALLOCATE( iproc1(inowin) ) 322 ALLOCATE( iobsi2(inowin) ) 323 ALLOCATE( iobsj2(inowin) ) 324 ALLOCATE( iproc2(inowin) ) 323 325 inowin = 0 324 326 DO ji = 1, inpfiles(jj)%nobs … … 334 336 END DO 335 337 336 CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' ) 338 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 339 CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 340 & iproc1, 'T' ) 341 iobsi2(:) = iobsi1(:) 342 iobsj2(:) = iobsj1(:) 343 iproc2(:) = iproc1(:) 344 ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 345 CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 346 & iproc1, 'U' ) 347 CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, & 348 & iproc2, 'V' ) 349 ENDIF 337 350 338 351 inowin = 0 … … 344 357 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 345 358 inowin = inowin + 1 346 inpfiles(jj)%iproc(ji,1) = iproc(inowin) 347 inpfiles(jj)%iobsi(ji,1) = iobsi(inowin) 348 inpfiles(jj)%iobsj(ji,1) = iobsj(inowin) 359 inpfiles(jj)%iproc(ji,1) = iproc1(inowin) 360 inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin) 361 inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin) 362 inpfiles(jj)%iproc(ji,2) = iproc2(inowin) 363 inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin) 364 inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin) 365 IF ( inpfiles(jj)%iproc(ji,1) /= & 366 & inpfiles(jj)%iproc(ji,2) ) THEN 367 CALL ctl_stop( 'Error in obs_read_prof:', & 368 & 'var1 and var2 observation on different processors') 369 ENDIF 349 370 ENDIF 350 371 END DO 351 DEALLOCATE( zlam, zphi, iobsi , iobsj, iproc)372 DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 ) 352 373 353 374 DO ji = 1, inpfiles(jj)%nobs … … 363 384 ENDIF 364 385 llvalprof = .FALSE. 365 IF ( ld t3d) THEN386 IF ( ldvar1 ) THEN 366 387 loop_t_count : DO ij = 1,inpfiles(jj)%nlev 367 388 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & … … 369 390 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 370 391 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 371 i t3dt0 = it3dt0 + 1392 ivar1t0 = ivar1t0 + 1 372 393 ENDIF 373 394 END DO loop_t_count 374 395 ENDIF 375 IF ( ld s3d) THEN396 IF ( ldvar2 ) THEN 376 397 loop_s_count : DO ij = 1,inpfiles(jj)%nlev 377 398 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & … … 379 400 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 380 401 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 381 i s3dt0 = is3dt0 + 1402 ivar2t0 = ivar2t0 + 1 382 403 ENDIF 383 404 END DO loop_s_count … … 388 409 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 389 410 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 390 & ld t3d) .OR. &411 & ldvar1 ) .OR. & 391 412 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 392 413 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 393 & ld s3d) ) THEN414 & ldvar2 ) ) THEN 394 415 ip3dt = ip3dt + 1 395 416 llvalprof = .TRUE. … … 405 426 406 427 END DO prof_files 407 428 408 429 !----------------------------------------------------------------------- 409 430 ! Get the time ordered indices of the input data … … 446 467 & zdat, & 447 468 & iindx ) 448 469 449 470 iv3dt(:) = -1 450 471 IF (ldsatt) THEN … … 452 473 iv3dt(2) = ip3dt 453 474 ELSE 454 iv3dt(1) = i t3dt0455 iv3dt(2) = i s3dt0475 iv3dt(1) = ivar1t0 476 iv3dt(2) = ivar2t0 456 477 ENDIF 457 478 CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 458 479 & kstp, jpi, jpj, jpk ) 459 480 460 481 ! * Read obs/positions, QC, all variable and assign to profdata 461 482 462 483 profdata%nprof = 0 463 484 profdata%nvprot(:) = 0 464 485 profdata%cvars(:) = clvars(:) 465 486 iprof = 0 466 487 467 488 ip3dt = 0 468 i t3dt = 0469 i s3dt = 0470 ityp t(:) = 0471 ityp tmpp(:) = 0472 473 ityp s(:) = 0474 ityp smpp(:) = 0475 476 ioserrcount = 0 489 ivar1t = 0 490 ivar2t = 0 491 itypvar1 (:) = 0 492 itypvar1mpp(:) = 0 493 494 itypvar2 (:) = 0 495 itypvar2mpp(:) = 0 496 497 ioserrcount = 0 477 498 DO jk = 1, iproftot 478 499 479 500 jj = ifileidx(iindx(jk)) 480 501 ji = iprofidx(iindx(jk)) … … 486 507 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 487 508 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 488 509 489 510 IF ( nproc == 0 ) THEN 490 511 IF ( inpfiles(jj)%iproc(ji,1) > nproc ) CYCLE … … 492 513 IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE 493 514 ENDIF 494 515 495 516 llvalprof = .FALSE. 496 517 … … 501 522 502 523 loop_prof : DO ij = 1, inpfiles(jj)%nlev 503 524 504 525 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 505 526 & CYCLE 506 527 507 528 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 508 529 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 509 530 510 531 llvalprof = .TRUE. 511 532 EXIT loop_prof 512 533 513 534 ENDIF 514 535 515 536 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 516 537 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 517 538 518 539 llvalprof = .TRUE. 519 540 EXIT loop_prof 520 541 521 542 ENDIF 522 543 523 544 END DO loop_prof 524 545 525 546 ! Set profile information 526 547 527 548 IF ( llvalprof ) THEN 528 549 529 550 iprof = iprof + 1 530 551 … … 545 566 profdata%nhou(iprof) = ihou 546 567 profdata%nmin(iprof) = imin 547 568 548 569 ! Profile space coordinates 549 570 profdata%rlam(iprof) = inpfiles(jj)%plam(ji) … … 551 572 552 573 ! Coordinate search parameters 553 profdata%mi (iprof,:) = inpfiles(jj)%iobsi(ji,1) 554 profdata%mj (iprof,:) = inpfiles(jj)%iobsj(ji,1) 555 574 profdata%mi (iprof,1) = inpfiles(jj)%iobsi(ji,1) 575 profdata%mj (iprof,1) = inpfiles(jj)%iobsj(ji,1) 576 profdata%mi (iprof,2) = inpfiles(jj)%iobsi(ji,2) 577 profdata%mj (iprof,2) = inpfiles(jj)%iobsj(ji,2) 578 556 579 ! Profile WMO number 557 580 profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji) 558 581 559 582 ! Instrument type 560 583 READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype … … 564 587 itype = 0 565 588 ENDIF 566 589 567 590 profdata%ntyp(iprof) = itype 568 591 569 592 ! QC stuff 570 593 … … 585 608 profdata%nqc(iprof) = 0 !TODO 586 609 587 loop_p : DO ij = 1, inpfiles(jj)%nlev 588 610 loop_p : DO ij = 1, inpfiles(jj)%nlev 611 589 612 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 590 613 & CYCLE … … 594 617 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 595 618 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 596 & ld t3d) .OR. &619 & ldvar1 ) .OR. & 597 620 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 598 621 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 599 & ld s3d) ) THEN622 & ldvar2 ) ) THEN 600 623 ip3dt = ip3dt + 1 601 624 ELSE 602 625 CYCLE 603 626 ENDIF 604 627 605 628 ENDIF 606 629 607 630 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 608 631 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 609 & ld t3d) .OR. ldsatt ) THEN610 632 & ldvar1 ) .OR. ldsatt ) THEN 633 611 634 IF (ldsatt) THEN 612 635 613 i t3dt = ip3dt636 ivar1t = ip3dt 614 637 615 638 ELSE 616 639 617 i t3dt = it3dt + 1618 640 ivar1t = ivar1t + 1 641 619 642 ENDIF 620 643 621 ! Depth of Tobservation622 profdata%var(1)%vdep(i t3dt) = &644 ! Depth of var1 observation 645 profdata%var(1)%vdep(ivar1t) = & 623 646 & inpfiles(jj)%pdep(ij,ji) 624 625 ! Depth of Tobservation QC626 profdata%var(1)%idqc(i t3dt) = &647 648 ! Depth of var1 observation QC 649 profdata%var(1)%idqc(ivar1t) = & 627 650 & inpfiles(jj)%idqc(ij,ji) 628 629 ! Depth of Tobservation QC flags630 profdata%var(1)%idqcf(:,i t3dt) = &651 652 ! Depth of var1 observation QC flags 653 profdata%var(1)%idqcf(:,ivar1t) = & 631 654 & inpfiles(jj)%idqcf(:,ij,ji) 632 655 633 656 ! Profile index 634 profdata%var(1)%nvpidx(i t3dt) = iprof635 657 profdata%var(1)%nvpidx(ivar1t) = iprof 658 636 659 ! Vertical index in original profile 637 profdata%var(1)%nvlidx(i t3dt) = ij638 639 ! Profile potential Tvalue660 profdata%var(1)%nvlidx(ivar1t) = ij 661 662 ! Profile var1 value 640 663 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 641 664 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 642 profdata%var(1)%vobs(i t3dt) = &665 profdata%var(1)%vobs(ivar1t) = & 643 666 & inpfiles(jj)%pob(ij,ji,1) 644 667 IF ( ldmod ) THEN 645 profdata%var(1)%vmod(i t3dt) = &668 profdata%var(1)%vmod(ivar1t) = & 646 669 & inpfiles(jj)%padd(ij,ji,1,1) 647 670 ENDIF 648 ! Count number of profile Tdata as function of type649 ityp t( profdata%ntyp(iprof) + 1 ) = &650 & ityp t( profdata%ntyp(iprof) + 1 ) + 1671 ! Count number of profile var1 data as function of type 672 itypvar1( profdata%ntyp(iprof) + 1 ) = & 673 & itypvar1( profdata%ntyp(iprof) + 1 ) + 1 651 674 ELSE 652 profdata%var(1)%vobs(i t3dt) = fbrmdi675 profdata%var(1)%vobs(ivar1t) = fbrmdi 653 676 ENDIF 654 677 655 ! Profile Tqc656 profdata%var(1)%nvqc(i t3dt) = &678 ! Profile var1 qc 679 profdata%var(1)%nvqc(ivar1t) = & 657 680 & inpfiles(jj)%ivlqc(ij,ji,1) 658 681 659 ! Profile Tqc flags660 profdata%var(1)%nvqcf(:,i t3dt) = &682 ! Profile var1 qc flags 683 profdata%var(1)%nvqcf(:,ivar1t) = & 661 684 & inpfiles(jj)%ivlqcf(:,ij,ji,1) 662 685 663 686 ! Profile insitu T value 664 profdata%var(1)%vext(it3dt,1) = & 665 & inpfiles(jj)%pext(ij,ji,1) 666 667 ENDIF 668 687 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 688 profdata%var(1)%vext(ivar1t,1) = & 689 & inpfiles(jj)%pext(ij,ji,1) 690 ENDIF 691 692 ENDIF 693 669 694 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 670 695 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 671 & ld s3d) .OR. ldsatt ) THEN672 696 & ldvar2 ) .OR. ldsatt ) THEN 697 673 698 IF (ldsatt) THEN 674 699 675 i s3dt = ip3dt700 ivar2t = ip3dt 676 701 677 702 ELSE 678 703 679 i s3dt = is3dt + 1680 704 ivar2t = ivar2t + 1 705 681 706 ENDIF 682 707 683 ! Depth of Sobservation684 profdata%var(2)%vdep(i s3dt) = &708 ! Depth of var2 observation 709 profdata%var(2)%vdep(ivar2t) = & 685 710 & inpfiles(jj)%pdep(ij,ji) 686 687 ! Depth of Sobservation QC688 profdata%var(2)%idqc(i s3dt) = &711 712 ! Depth of var2 observation QC 713 profdata%var(2)%idqc(ivar2t) = & 689 714 & inpfiles(jj)%idqc(ij,ji) 690 691 ! Depth of Sobservation QC flags692 profdata%var(2)%idqcf(:,i s3dt) = &715 716 ! Depth of var2 observation QC flags 717 profdata%var(2)%idqcf(:,ivar2t) = & 693 718 & inpfiles(jj)%idqcf(:,ij,ji) 694 719 695 720 ! Profile index 696 profdata%var(2)%nvpidx(i s3dt) = iprof697 721 profdata%var(2)%nvpidx(ivar2t) = iprof 722 698 723 ! Vertical index in original profile 699 profdata%var(2)%nvlidx(i s3dt) = ij700 701 ! Profile Svalue724 profdata%var(2)%nvlidx(ivar2t) = ij 725 726 ! Profile var2 value 702 727 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 703 728 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 704 profdata%var(2)%vobs(i s3dt) = &729 profdata%var(2)%vobs(ivar2t) = & 705 730 & inpfiles(jj)%pob(ij,ji,2) 706 731 IF ( ldmod ) THEN 707 profdata%var(2)%vmod(i s3dt) = &732 profdata%var(2)%vmod(ivar2t) = & 708 733 & inpfiles(jj)%padd(ij,ji,1,2) 709 734 ENDIF 710 ! Count number of profile Sdata as function of type711 ityp s( profdata%ntyp(iprof) + 1 ) = &712 & ityp s( profdata%ntyp(iprof) + 1 ) + 1735 ! Count number of profile var2 data as function of type 736 itypvar2( profdata%ntyp(iprof) + 1 ) = & 737 & itypvar2( profdata%ntyp(iprof) + 1 ) + 1 713 738 ELSE 714 profdata%var(2)%vobs(i s3dt) = fbrmdi739 profdata%var(2)%vobs(ivar2t) = fbrmdi 715 740 ENDIF 716 717 ! Profile Sqc718 profdata%var(2)%nvqc(i s3dt) = &741 742 ! Profile var2 qc 743 profdata%var(2)%nvqc(ivar2t) = & 719 744 & inpfiles(jj)%ivlqc(ij,ji,2) 720 745 721 ! Profile Sqc flags722 profdata%var(2)%nvqcf(:,i s3dt) = &746 ! Profile var2 qc flags 747 profdata%var(2)%nvqcf(:,ivar2t) = & 723 748 & inpfiles(jj)%ivlqcf(:,ij,ji,2) 724 749 725 750 ENDIF 726 751 727 752 END DO loop_p 728 753 … … 736 761 ! Sum up over processors 737 762 !----------------------------------------------------------------------- 738 739 CALL obs_mpp_sum_integer ( i t3dt0, it3dtmpp )740 CALL obs_mpp_sum_integer ( i s3dt0, is3dtmpp )741 CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp)742 743 CALL obs_mpp_sum_integers( ityp t, ityptmpp, ntyp1770 + 1 )744 CALL obs_mpp_sum_integers( ityp s, itypsmpp, ntyp1770 + 1 )745 763 764 CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp ) 765 CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 766 CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp ) 767 768 CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 ) 769 CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 770 746 771 !----------------------------------------------------------------------- 747 772 ! Output number of observations. … … 749 774 IF(lwp) THEN 750 775 WRITE(numout,*) 751 WRITE(numout,'( 1X,A)') 'Profile data'776 WRITE(numout,'(A)') ' Profile data' 752 777 WRITE(numout,'(1X,A)') '------------' 753 778 WRITE(numout,*) 754 WRITE(numout,'(1X,A)') 'Profile T data'755 WRITE(numout,'(1X,A)') '-------------- '779 WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) ) 780 WRITE(numout,'(1X,A)') '------------------------' 756 781 DO ji = 0, ntyp1770 757 IF ( ityp tmpp(ji+1) > 0 ) THEN782 IF ( itypvar1mpp(ji+1) > 0 ) THEN 758 783 WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 759 784 & cwmonam1770(ji)(1:52),' = ', & 760 & ityp tmpp(ji+1)785 & itypvar1mpp(ji+1) 761 786 ENDIF 762 787 END DO … … 764 789 & '---------------------------------------------------------------' 765 790 WRITE(numout,'(1X,A55,I8)') & 766 & 'Total profile T data = ',&767 & it3dtmpp791 & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// & 792 & ' = ', ivar1tmpp 768 793 WRITE(numout,'(1X,A)') & 769 794 & '---------------------------------------------------------------' 770 795 WRITE(numout,*) 771 WRITE(numout,'(1X,A)') 'Profile S data'772 WRITE(numout,'(1X,A)') '-------------- '796 WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) ) 797 WRITE(numout,'(1X,A)') '------------------------' 773 798 DO ji = 0, ntyp1770 774 IF ( ityp smpp(ji+1) > 0 ) THEN799 IF ( itypvar2mpp(ji+1) > 0 ) THEN 775 800 WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 776 801 & cwmonam1770(ji)(1:52),' = ', & 777 & ityp smpp(ji+1)802 & itypvar2mpp(ji+1) 778 803 ENDIF 779 804 END DO … … 781 806 & '---------------------------------------------------------------' 782 807 WRITE(numout,'(1X,A55,I8)') & 783 & 'Total profile S data = ',&784 & is3dtmpp808 & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// & 809 & ' = ', ivar2tmpp 785 810 WRITE(numout,'(1X,A)') & 786 811 & '---------------------------------------------------------------' 787 812 WRITE(numout,*) 788 813 ENDIF 789 814 790 815 IF (ldsatt) THEN 791 816 profdata%nvprot(1) = ip3dt … … 794 819 profdata%nvprotmpp(2) = ip3dtmpp 795 820 ELSE 796 profdata%nvprot(1) = i t3dt797 profdata%nvprot(2) = i s3dt798 profdata%nvprotmpp(1) = i t3dtmpp799 profdata%nvprotmpp(2) = i s3dtmpp821 profdata%nvprot(1) = ivar1t 822 profdata%nvprot(2) = ivar2t 823 profdata%nvprotmpp(1) = ivar1tmpp 824 profdata%nvprotmpp(2) = ivar2tmpp 800 825 ENDIF 801 826 profdata%nprof = iprof … … 804 829 ! Model level search 805 830 !----------------------------------------------------------------------- 806 IF ( ld t3d) THEN831 IF ( ldvar1 ) THEN 807 832 CALL obs_level_search( jpk, gdept_1d, & 808 833 & profdata%nvprot(1), profdata%var(1)%vdep, & 809 834 & profdata%var(1)%mvk ) 810 835 ENDIF 811 IF ( ld s3d) THEN836 IF ( ldvar2 ) THEN 812 837 CALL obs_level_search( jpk, gdept_1d, & 813 838 & profdata%nvprot(2), profdata%var(2)%vdep, & 814 839 & profdata%var(2)%mvk ) 815 840 ENDIF 816 841 817 842 !----------------------------------------------------------------------- 818 843 ! Set model equivalent to 99999 … … 826 851 ! Deallocate temporary data 827 852 !----------------------------------------------------------------------- 828 DEALLOCATE( ifileidx, iprofidx, zdat )853 DEALLOCATE( ifileidx, iprofidx, zdat, clvars ) 829 854 830 855 !----------------------------------------------------------------------- … … 836 861 DEALLOCATE( inpfiles ) 837 862 838 END SUBROUTINE obs_rea_pro _dri863 END SUBROUTINE obs_rea_prof 839 864 840 865 END MODULE obs_read_prof
Note: See TracChangeset
for help on using the changeset viewer.