New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r4292 r6225  
    2525   USE netcdf                   ! NetCDF library 
    2626   USE obs_oper                 ! Observation operators 
    27    USE obs_prof_io              ! Profile files I/O (non-FB files) 
    2827   USE lib_mpp                  ! For ctl_warn/stop 
     28   USE obs_fbm                  ! Feedback routines 
    2929 
    3030   IMPLICIT NONE 
    … …  
    3333   PRIVATE 
    3434 
    35    PUBLIC obs_rea_pro_dri  ! Read the profile observations  
     35   PUBLIC obs_rea_prof  ! Read the profile observations  
    3636 
    3737   !!---------------------------------------------------------------------- 
    … …  
    4242 
    4343CONTAINS 
    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 ) 
    5049      !!--------------------------------------------------------------------- 
    5150      !! 
    52       !!                   *** ROUTINE obs_rea_pro_dri *** 
     51      !!                   *** ROUTINE obs_rea_prof *** 
    5352      !! 
    5453      !! ** Purpose : Read from file the profile observations 
    5554      !! 
    56       !! ** Method  : Depending on kformat either ENACT, CORIOLIS or 
    57       !!              feedback data files are read 
     55      !! ** Method  : Read feedback data in and transform to NEMO internal  
     56      !!              profile data structure 
    5857      !! 
    5958      !! ** Action  :  
    … …  
    6362      !! History :   
    6463      !!      ! :  2009-09 (K. Mogensen) : New merged version of old routines 
     64      !!      ! :  2015-08 (M. Martin) : Merged profile and velocity routines 
    6565      !!---------------------------------------------------------------------- 
    66       !! * Modules used 
    67     
     66 
    6867      !! * 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 
    7471      CHARACTER(LEN=128), INTENT(IN) ::  & 
    75          & cfilenames(knumfiles)  ! File names to read in 
     72         & cdfilenames(knumfiles)        ! File names to read in 
    7673      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 
    8783      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    88          & kdailyavtypes 
     84         & kdailyavtypes                ! Types of daily average observations 
    8985 
    9086      !! * 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 
    9290      INTEGER :: jvar 
    9391      INTEGER :: ji 
    … …  
    105103      INTEGER :: imin 
    106104      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 
    107118      INTEGER, DIMENSION(knumfiles) :: & 
    108119         & irefdate 
    109120      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  
    117125      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
    118          & iobsi,    & 
    119          & iobsj,    & 
    120          & iproc,    & 
     126         & iobsi1,    & 
     127         & iobsj1,    & 
     128         & iproc1,    & 
     129         & iobsi2,    & 
     130         & iobsj2,    & 
     131         & iproc2,    & 
    121132         & iindx,    & 
    122133         & ifileidx, & 
    123134         & iprofidx 
    124       INTEGER :: itype 
    125135      INTEGER, DIMENSION(imaxavtypes) :: & 
    126136         & idailyavtypes 
     137      INTEGER, DIMENSION(kvars) :: & 
     138         & iv3dt 
    127139      REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
    128140         & zphi, & 
    129141         & zlam 
    130       REAL(dp), DIMENSION(:), ALLOCATABLE :: & 
     142      REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
    131143         & zdat 
     144      REAL(wp), DIMENSION(knumfiles) :: & 
     145         & djulini, & 
     146         & djulend 
    132147      LOGICAL :: llvalprof 
     148      LOGICAL :: lldavtimset 
    133149      TYPE(obfbdata), POINTER, DIMENSION(:) :: & 
    134150         & inpfiles 
    135       REAL(dp), 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, DIMENSION(kvars) :: & 
    146          & iv3dt 
    147       CHARACTER(len=8) :: cl_refdate 
    148     
     151 
    149152      ! Local initialization 
    150153      iprof = 0 
    151       it3dt0 = 0 
    152       is3dt0 = 0 
     154      ivar1t0 = 0 
     155      ivar2t0 = 0 
    153156      ip3dt = 0 
    154157 
    155158      ! Daily average types 
     159      lldavtimset = .FALSE. 
    156160      IF ( PRESENT(kdailyavtypes) ) THEN 
    157161         idailyavtypes(:) = kdailyavtypes(:) 
     162         IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE. 
    158163      ELSE 
    159164         idailyavtypes(:) = -1 
    … …  
    161166 
    162167      !----------------------------------------------------------------------- 
    163       ! Check data the model part is just with feedback data files 
    164       !----------------------------------------------------------------------- 
    165       IF ( ldmod .AND. ( kformat /= 0 ) ) THEN 
    166          CALL ctl_stop( 'Model can only be read from feedback data' ) 
    167          RETURN 
    168       ENDIF 
    169  
    170       !----------------------------------------------------------------------- 
    171168      ! Count the number of files needed and allocate the obfbdata type 
    172169      !----------------------------------------------------------------------- 
    173        
     170 
    174171      inobf = knumfiles 
    175        
     172 
    176173      ALLOCATE( inpfiles(inobf) ) 
    177174 
    178175      prof_files : DO jj = 1, inobf 
    179            
     176 
    180177         !--------------------------------------------------------------------- 
    181178         ! Prints 
    … …  
    184181            WRITE(numout,*) 
    185182            WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', & 
    186                & TRIM( TRIM( cfilenames(jj) ) ) 
     183               & TRIM( TRIM( cdfilenames(jj) ) ) 
    187184            WRITE(numout,*) ' ~~~~~~~~~~~~~~~' 
    188185            WRITE(numout,*) 
    … …  
    192189         !  Initialization: Open file and get dimensions only 
    193190         !--------------------------------------------------------------------- 
    194           
    195          iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, & 
     191 
     192         iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, & 
    196193            &                      i_file_id ) 
    197           
     194 
    198195         IF ( iflag /= nf90_noerr ) THEN 
    199196 
    200197            IF ( ldignmis ) THEN 
    201198               inpfiles(jj)%nobs = 0 
    202                CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // & 
     199               CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // & 
    203200                  &           ' not found' ) 
    204201            ELSE  
    205                CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // & 
     202               CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // & 
    206203                  &           ' not found' ) 
    207204            ENDIF 
    208205 
    209206         ELSE  
    210              
     207 
    211208            !------------------------------------------------------------------ 
    212             !  Close the file since it is opened in read_proffile 
     209            !  Close the file since it is opened in read_obfbdata 
    213210            !------------------------------------------------------------------ 
    214              
     211 
    215212            iflag = nf90_close( i_file_id ) 
    216213 
    … …  
    218215            !  Read the profile file into inpfiles 
    219216            !------------------------------------------------------------------ 
    220             IF ( kformat == 0 ) THEN 
    221                CALL init_obfbdata( inpfiles(jj) ) 
    222                IF(lwp) THEN 
    223                   WRITE(numout,*) 
    224                   WRITE(numout,*)'Reading from feedback file :', & 
    225                      &           TRIM( cfilenames(jj) ) 
    226                ENDIF 
    227                CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), & 
    228                   &                ldgrid = .TRUE. ) 
    229                IF ( inpfiles(jj)%nvar < 2 ) THEN 
    230                   CALL ctl_stop( 'Feedback format error' ) 
    231                   RETURN 
    232                ENDIF 
    233                IF ( TRIM(inpfiles(jj)%cname(1)) /= 'POTM' ) THEN 
    234                   CALL ctl_stop( 'Feedback format error' ) 
    235                   RETURN 
    236                ENDIF 
    237                IF ( TRIM(inpfiles(jj)%cname(2)) /= 'PSAL' ) THEN 
    238                   CALL ctl_stop( 'Feedback format error' ) 
    239                   RETURN 
    240                ENDIF 
    241                IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 
    242                   CALL ctl_stop( 'Model not in input data' ) 
    243                   RETURN 
    244                ENDIF 
    245             ELSEIF ( kformat == 1 ) THEN 
    246                CALL read_enactfile( TRIM( cfilenames(jj) ), inpfiles(jj), & 
    247                   &                 numout, lwp, .TRUE. ) 
    248             ELSEIF ( kformat == 2 ) THEN 
    249                CALL read_coriofile( TRIM( cfilenames(jj) ), inpfiles(jj), & 
    250                   &                 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 
    251235            ELSE 
    252                CALL ctl_stop( 'File format unknown' ) 
    253             ENDIF 
    254              
     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 
    255244            !------------------------------------------------------------------ 
    256245            !  Change longitude (-180,180) 
    … …  
    270259            !  Calculate the date  (change eventually) 
    271260            !------------------------------------------------------------------ 
    272             cl_refdate=inpfiles(jj)%cdjuldref(1:8) 
    273             READ(cl_refdate,'(I8)') irefdate(jj) 
    274              
     261            clrefdate=inpfiles(jj)%cdjuldref(1:8) 
     262            READ(clrefdate,'(I8)') irefdate(jj) 
     263 
    275264            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec ) 
    276265            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), & 
    … …  
    280269               &           krefdate = irefdate(jj) ) 
    281270 
    282             IF ( ldavtimset ) THEN 
     271            ioserrcount=0 
     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 
    283279               DO ji = 1, inpfiles(jj)%nobs 
    284                   !  
    285                   !  for daily averaged data for example 
    286                   !  MRB data (itype==820) force the time 
    287                   !  to be the  end of the day 
    288                   ! 
    289                   READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype 
    290                   IF ( ANY (idailyavtypes == itype ) ) THEN 
    291                      inpfiles(jj)%ptim(ji) = & 
    292                      & INT(inpfiles(jj)%ptim(ji)) + 1 
    293                   ENDIF 
     280                  READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype 
     281900               IF ( ios /= 0 ) THEN 
     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 
    294298               END DO 
    295             ENDIF 
    296              
     299 
     300            ENDIF 
     301 
    297302            IF ( inpfiles(jj)%nobs > 0 ) THEN 
    298                inpfiles(jj)%iproc = -1 
    299                inpfiles(jj)%iobsi = -1 
    300                inpfiles(jj)%iobsj = -1 
     303               inpfiles(jj)%iproc(:,:) = -1 
     304               inpfiles(jj)%iobsi(:,:) = -1 
     305               inpfiles(jj)%iobsj(:,:) = -1 
    301306            ENDIF 
    302307            inowin = 0 
    … …  
    312317            ALLOCATE( zlam(inowin)  ) 
    313318            ALLOCATE( zphi(inowin)  ) 
    314             ALLOCATE( iobsi(inowin) ) 
    315             ALLOCATE( iobsj(inowin) ) 
    316             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) ) 
    317325            inowin = 0 
    318326            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    328336            END DO 
    329337 
    330             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 
    331350 
    332351            inowin = 0 
    … …  
    338357                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    339358                  inowin = inowin + 1 
    340                   inpfiles(jj)%iproc(ji,1) = iproc(inowin) 
    341                   inpfiles(jj)%iobsi(ji,1) = iobsi(inowin) 
    342                   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 
    343370               ENDIF 
    344371            END DO 
    345             DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc ) 
     372            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 ) 
    346373 
    347374            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    357384                  ENDIF 
    358385                  llvalprof = .FALSE. 
    359                   IF ( ldt3d ) THEN 
     386                  IF ( ldvar1 ) THEN 
    360387                     loop_t_count : DO ij = 1,inpfiles(jj)%nlev 
    361388                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    … …  
    363390                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    364391                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    365                            it3dt0 = it3dt0 + 1 
     392                           ivar1t0 = ivar1t0 + 1 
    366393                        ENDIF 
    367394                     END DO loop_t_count 
    368395                  ENDIF 
    369                   IF ( lds3d ) THEN 
     396                  IF ( ldvar2 ) THEN 
    370397                     loop_s_count : DO ij = 1,inpfiles(jj)%nlev 
    371398                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    … …  
    373400                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    374401                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    375                            is3dt0 = is3dt0 + 1 
     402                           ivar2t0 = ivar2t0 + 1 
    376403                        ENDIF 
    377404                     END DO loop_s_count 
    … …  
    382409                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    383410                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    384                         &     ldt3d ) .OR. & 
     411                        &     ldvar1 ) .OR. & 
    385412                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    386413                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    387                         &     lds3d ) ) THEN 
     414                        &     ldvar2 ) ) THEN 
    388415                        ip3dt = ip3dt + 1 
    389416                        llvalprof = .TRUE. 
    … …  
    399426 
    400427      END DO prof_files 
    401        
     428 
    402429      !----------------------------------------------------------------------- 
    403430      ! Get the time ordered indices of the input data 
    … …  
    440467         &               zdat,     & 
    441468         &               iindx   ) 
    442        
     469 
    443470      iv3dt(:) = -1 
    444471      IF (ldsatt) THEN 
    … …  
    446473         iv3dt(2) = ip3dt 
    447474      ELSE 
    448          iv3dt(1) = it3dt0 
    449          iv3dt(2) = is3dt0 
     475         iv3dt(1) = ivar1t0 
     476         iv3dt(2) = ivar2t0 
    450477      ENDIF 
    451478      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 
    452479         &                 kstp, jpi, jpj, jpk ) 
    453        
     480 
    454481      ! * Read obs/positions, QC, all variable and assign to profdata 
    455482 
    456483      profdata%nprof     = 0 
    457484      profdata%nvprot(:) = 0 
    458  
     485      profdata%cvars(:)  = clvars(:) 
    459486      iprof = 0 
    460487 
    461488      ip3dt = 0 
    462       it3dt = 0 
    463       is3dt = 0 
    464       itypt   (:) = 0 
    465       ityptmpp(:) = 0 
    466        
    467       ityps   (:) = 0 
    468       itypsmpp(:) = 0 
    469        
    470        
     489      ivar1t = 0 
     490      ivar2t = 0 
     491      itypvar1   (:) = 0 
     492      itypvar1mpp(:) = 0 
     493 
     494      itypvar2   (:) = 0 
     495      itypvar2mpp(:) = 0 
     496 
     497      ioserrcount = 0 
    471498      DO jk = 1, iproftot 
    472           
     499 
    473500         jj = ifileidx(iindx(jk)) 
    474501         ji = iprofidx(iindx(jk)) 
    … …  
    480507         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
    481508            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 
    482              
     509 
    483510            IF ( nproc == 0 ) THEN 
    484511               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE 
    … …  
    486513               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE 
    487514            ENDIF 
    488              
     515 
    489516            llvalprof = .FALSE. 
    490517 
    … …  
    495522 
    496523            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
    497                 
     524 
    498525               IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    499526                  & CYCLE 
    500                 
     527 
    501528               IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    502529                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    503                    
     530 
    504531                  llvalprof = .TRUE.  
    505532                  EXIT loop_prof 
    506                    
     533 
    507534               ENDIF 
    508                 
     535 
    509536               IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    510537                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    511                    
     538 
    512539                  llvalprof = .TRUE.  
    513540                  EXIT loop_prof 
    514                    
     541 
    515542               ENDIF 
    516                 
     543 
    517544            END DO loop_prof 
    518              
     545 
    519546            ! Set profile information 
    520              
     547 
    521548            IF ( llvalprof ) THEN 
    522                 
     549 
    523550               iprof = iprof + 1 
    524551 
    … …  
    539566               profdata%nhou(iprof) = ihou 
    540567               profdata%nmin(iprof) = imin 
    541                 
     568 
    542569               ! Profile space coordinates 
    543570               profdata%rlam(iprof) = inpfiles(jj)%plam(ji) 
    … …  
    545572 
    546573               ! Coordinate search parameters 
    547                profdata%mi  (iprof,:) = inpfiles(jj)%iobsi(ji,1) 
    548                profdata%mj  (iprof,:) = inpfiles(jj)%iobsj(ji,1) 
    549                 
     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 
    550579               ! Profile WMO number 
    551580               profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji) 
    552                 
     581 
    553582               ! Instrument type 
    554                READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype 
     583               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype 
     584901            IF ( ios /= 0 ) THEN 
     585                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
     586                  ioserrcount = ioserrcount + 1 
     587                  itype = 0 
     588               ENDIF 
     589 
    555590               profdata%ntyp(iprof) = itype 
    556                 
     591 
    557592               ! QC stuff 
    558593 
    … …  
    573608               profdata%nqc(iprof)  = 0 !TODO 
    574609 
    575                loop_p : DO ij = 1, inpfiles(jj)%nlev             
    576                    
     610               loop_p : DO ij = 1, inpfiles(jj)%nlev 
     611 
    577612                  IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    578613                     & CYCLE 
    … …  
    582617                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    583618                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    584                         &     ldt3d ) .OR. & 
     619                        &     ldvar1 ) .OR. & 
    585620                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    586621                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    587                         &     lds3d ) ) THEN 
     622                        &     ldvar2 ) ) THEN 
    588623                        ip3dt = ip3dt + 1 
    589624                     ELSE 
    590625                        CYCLE 
    591626                     ENDIF 
    592                       
     627 
    593628                  ENDIF 
    594629 
    595630                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    596631                     &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    597                      &       ldt3d ) .OR. ldsatt ) THEN 
    598                       
     632                     &       ldvar1 ) .OR. ldsatt ) THEN 
     633 
    599634                     IF (ldsatt) THEN 
    600635 
    601                         it3dt = ip3dt 
     636                        ivar1t = ip3dt 
    602637 
    603638                     ELSE 
    604639 
    605                         it3dt = it3dt + 1 
    606                          
     640                        ivar1t = ivar1t + 1 
     641 
    607642                     ENDIF 
    608643 
    609                      ! Depth of T observation 
    610                      profdata%var(1)%vdep(it3dt) = & 
     644                     ! Depth of var1 observation 
     645                     profdata%var(1)%vdep(ivar1t) = & 
    611646                        &                inpfiles(jj)%pdep(ij,ji) 
    612                       
    613                      ! Depth of T observation QC 
    614                      profdata%var(1)%idqc(it3dt) = & 
     647 
     648                     ! Depth of var1 observation QC 
     649                     profdata%var(1)%idqc(ivar1t) = & 
    615650                        &                inpfiles(jj)%idqc(ij,ji) 
    616                       
    617                      ! Depth of T observation QC flags 
    618                      profdata%var(1)%idqcf(:,it3dt) = & 
     651 
     652                     ! Depth of var1 observation QC flags 
     653                     profdata%var(1)%idqcf(:,ivar1t) = & 
    619654                        &                inpfiles(jj)%idqcf(:,ij,ji) 
    620                       
     655 
    621656                     ! Profile index 
    622                      profdata%var(1)%nvpidx(it3dt) = iprof 
    623                       
     657                     profdata%var(1)%nvpidx(ivar1t) = iprof 
     658 
    624659                     ! Vertical index in original profile 
    625                      profdata%var(1)%nvlidx(it3dt) = ij 
    626  
    627                      ! Profile potential T value 
     660                     profdata%var(1)%nvlidx(ivar1t) = ij 
     661 
     662                     ! Profile var1 value 
    628663                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    629664                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    630                         profdata%var(1)%vobs(it3dt) = & 
     665                        profdata%var(1)%vobs(ivar1t) = & 
    631666                           &                inpfiles(jj)%pob(ij,ji,1) 
    632667                        IF ( ldmod ) THEN 
    633                            profdata%var(1)%vmod(it3dt) = & 
     668                           profdata%var(1)%vmod(ivar1t) = & 
    634669                              &                inpfiles(jj)%padd(ij,ji,1,1) 
    635670                        ENDIF 
    636                         ! Count number of profile T data as function of type 
    637                         itypt( profdata%ntyp(iprof) + 1 ) = & 
    638                            & itypt( profdata%ntyp(iprof) + 1 ) + 1 
     671                        ! Count number of profile var1 data as function of type 
     672                        itypvar1( profdata%ntyp(iprof) + 1 ) = & 
     673                           & itypvar1( profdata%ntyp(iprof) + 1 ) + 1 
    639674                     ELSE 
    640                         profdata%var(1)%vobs(it3dt) = fbrmdi 
     675                        profdata%var(1)%vobs(ivar1t) = fbrmdi 
    641676                     ENDIF 
    642677 
    643                      ! Profile T qc 
    644                      profdata%var(1)%nvqc(it3dt) = & 
     678                     ! Profile var1 qc 
     679                     profdata%var(1)%nvqc(ivar1t) = & 
    645680                        & inpfiles(jj)%ivlqc(ij,ji,1) 
    646681 
    647                      ! Profile T qc flags 
    648                      profdata%var(1)%nvqcf(:,it3dt) = & 
     682                     ! Profile var1 qc flags 
     683                     profdata%var(1)%nvqcf(:,ivar1t) = & 
    649684                        & inpfiles(jj)%ivlqcf(:,ij,ji,1) 
    650685 
    651686                     ! Profile insitu T value 
    652                      profdata%var(1)%vext(it3dt,1) = & 
    653                         &                inpfiles(jj)%pext(ij,ji,1) 
    654                       
    655                   ENDIF 
    656                    
     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 
    657694                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    658695                     &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    659                      &   lds3d ) .OR. ldsatt ) THEN 
    660                       
     696                     &   ldvar2 ) .OR. ldsatt ) THEN 
     697 
    661698                     IF (ldsatt) THEN 
    662699 
    663                         is3dt = ip3dt 
     700                        ivar2t = ip3dt 
    664701 
    665702                     ELSE 
    666703 
    667                         is3dt = is3dt + 1 
    668                          
     704                        ivar2t = ivar2t + 1 
     705 
    669706                     ENDIF 
    670707 
    671                      ! Depth of S observation 
    672                      profdata%var(2)%vdep(is3dt) = & 
     708                     ! Depth of var2 observation 
     709                     profdata%var(2)%vdep(ivar2t) = & 
    673710                        &                inpfiles(jj)%pdep(ij,ji) 
    674                       
    675                      ! Depth of S observation QC 
    676                      profdata%var(2)%idqc(is3dt) = & 
     711 
     712                     ! Depth of var2 observation QC 
     713                     profdata%var(2)%idqc(ivar2t) = & 
    677714                        &                inpfiles(jj)%idqc(ij,ji) 
    678                       
    679                      ! Depth of S observation QC flags 
    680                      profdata%var(2)%idqcf(:,is3dt) = & 
     715 
     716                     ! Depth of var2 observation QC flags 
     717                     profdata%var(2)%idqcf(:,ivar2t) = & 
    681718                        &                inpfiles(jj)%idqcf(:,ij,ji) 
    682                       
     719 
    683720                     ! Profile index 
    684                      profdata%var(2)%nvpidx(is3dt) = iprof 
    685                       
     721                     profdata%var(2)%nvpidx(ivar2t) = iprof 
     722 
    686723                     ! Vertical index in original profile 
    687                      profdata%var(2)%nvlidx(is3dt) = ij 
    688  
    689                      ! Profile S value 
     724                     profdata%var(2)%nvlidx(ivar2t) = ij 
     725 
     726                     ! Profile var2 value 
    690727                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    691728                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    692                         profdata%var(2)%vobs(is3dt) = & 
     729                        profdata%var(2)%vobs(ivar2t) = & 
    693730                           &                inpfiles(jj)%pob(ij,ji,2) 
    694731                        IF ( ldmod ) THEN 
    695                            profdata%var(2)%vmod(is3dt) = & 
     732                           profdata%var(2)%vmod(ivar2t) = & 
    696733                              &                inpfiles(jj)%padd(ij,ji,1,2) 
    697734                        ENDIF 
    698                         ! Count number of profile S data as function of type 
    699                         ityps( profdata%ntyp(iprof) + 1 ) = & 
    700                            & ityps( profdata%ntyp(iprof) + 1 ) + 1 
     735                        ! Count number of profile var2 data as function of type 
     736                        itypvar2( profdata%ntyp(iprof) + 1 ) = & 
     737                           & itypvar2( profdata%ntyp(iprof) + 1 ) + 1 
    701738                     ELSE 
    702                         profdata%var(2)%vobs(is3dt) = fbrmdi 
     739                        profdata%var(2)%vobs(ivar2t) = fbrmdi 
    703740                     ENDIF 
    704                       
    705                      ! Profile S qc 
    706                      profdata%var(2)%nvqc(is3dt) = & 
     741 
     742                     ! Profile var2 qc 
     743                     profdata%var(2)%nvqc(ivar2t) = & 
    707744                        & inpfiles(jj)%ivlqc(ij,ji,2) 
    708745 
    709                      ! Profile S qc flags 
    710                      profdata%var(2)%nvqcf(:,is3dt) = & 
     746                     ! Profile var2 qc flags 
     747                     profdata%var(2)%nvqcf(:,ivar2t) = & 
    711748                        & inpfiles(jj)%ivlqcf(:,ij,ji,2) 
    712749 
    713750                  ENDIF 
    714              
     751 
    715752               END DO loop_p 
    716753 
    … …  
    724761      ! Sum up over processors 
    725762      !----------------------------------------------------------------------- 
    726        
    727       CALL obs_mpp_sum_integer ( it3dt0, it3dtmpp ) 
    728       CALL obs_mpp_sum_integer ( is3dt0, is3dtmpp ) 
    729       CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp ) 
    730        
    731       CALL obs_mpp_sum_integers( itypt, ityptmpp, ntyp1770 + 1 ) 
    732       CALL obs_mpp_sum_integers( ityps, itypsmpp, ntyp1770 + 1 ) 
    733        
     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 
    734771      !----------------------------------------------------------------------- 
    735772      ! Output number of observations. 
    … …  
    737774      IF(lwp) THEN 
    738775         WRITE(numout,*)  
    739          WRITE(numout,'(1X,A)') 'Profile data' 
     776         WRITE(numout,'(A)') ' Profile data' 
    740777         WRITE(numout,'(1X,A)') '------------' 
    741778         WRITE(numout,*)  
    742          WRITE(numout,'(1X,A)') 'Profile T data' 
    743          WRITE(numout,'(1X,A)') '--------------' 
     779         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) ) 
     780         WRITE(numout,'(1X,A)') '------------------------' 
    744781         DO ji = 0, ntyp1770 
    745             IF ( ityptmpp(ji+1) > 0 ) THEN 
     782            IF ( itypvar1mpp(ji+1) > 0 ) THEN 
    746783               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
    747784                  & cwmonam1770(ji)(1:52),' = ', & 
    748                   & ityptmpp(ji+1) 
     785                  & itypvar1mpp(ji+1) 
    749786            ENDIF 
    750787         END DO 
    … …  
    752789            & '---------------------------------------------------------------' 
    753790         WRITE(numout,'(1X,A55,I8)') & 
    754             & 'Total profile T data                                 = ',& 
    755             & it3dtmpp 
     791            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// & 
     792            & '             = ', ivar1tmpp 
    756793         WRITE(numout,'(1X,A)') & 
    757794            & '---------------------------------------------------------------' 
    758795         WRITE(numout,*)  
    759          WRITE(numout,'(1X,A)') 'Profile S data' 
    760          WRITE(numout,'(1X,A)') '--------------' 
     796         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) ) 
     797         WRITE(numout,'(1X,A)') '------------------------' 
    761798         DO ji = 0, ntyp1770 
    762             IF ( itypsmpp(ji+1) > 0 ) THEN 
     799            IF ( itypvar2mpp(ji+1) > 0 ) THEN 
    763800               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
    764801                  & cwmonam1770(ji)(1:52),' = ', & 
    765                   & itypsmpp(ji+1) 
     802                  & itypvar2mpp(ji+1) 
    766803            ENDIF 
    767804         END DO 
    … …  
    769806            & '---------------------------------------------------------------' 
    770807         WRITE(numout,'(1X,A55,I8)') & 
    771             & 'Total profile S data                                 = ',& 
    772             & is3dtmpp 
     808            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// & 
     809            & '             = ', ivar2tmpp 
    773810         WRITE(numout,'(1X,A)') & 
    774811            & '---------------------------------------------------------------' 
    775812         WRITE(numout,*)  
    776813      ENDIF 
    777        
     814 
    778815      IF (ldsatt) THEN 
    779816         profdata%nvprot(1)    = ip3dt 
    … …  
    782819         profdata%nvprotmpp(2) = ip3dtmpp 
    783820      ELSE 
    784          profdata%nvprot(1)    = it3dt 
    785          profdata%nvprot(2)    = is3dt 
    786          profdata%nvprotmpp(1) = it3dtmpp 
    787          profdata%nvprotmpp(2) = is3dtmpp 
     821         profdata%nvprot(1)    = ivar1t 
     822         profdata%nvprot(2)    = ivar2t 
     823         profdata%nvprotmpp(1) = ivar1tmpp 
     824         profdata%nvprotmpp(2) = ivar2tmpp 
    788825      ENDIF 
    789826      profdata%nprof        = iprof 
    … …  
    792829      ! Model level search 
    793830      !----------------------------------------------------------------------- 
    794       IF ( ldt3d ) THEN 
     831      IF ( ldvar1 ) THEN 
    795832         CALL obs_level_search( jpk, gdept_1d, & 
    796833            & profdata%nvprot(1), profdata%var(1)%vdep, & 
    797834            & profdata%var(1)%mvk ) 
    798835      ENDIF 
    799       IF ( lds3d ) THEN 
     836      IF ( ldvar2 ) THEN 
    800837         CALL obs_level_search( jpk, gdept_1d, & 
    801838            & profdata%nvprot(2), profdata%var(2)%vdep, & 
    802839            & profdata%var(2)%mvk ) 
    803840      ENDIF 
    804        
     841 
    805842      !----------------------------------------------------------------------- 
    806843      ! Set model equivalent to 99999 
    … …  
    814851      ! Deallocate temporary data 
    815852      !----------------------------------------------------------------------- 
    816       DEALLOCATE( ifileidx, iprofidx, zdat ) 
     853      DEALLOCATE( ifileidx, iprofidx, zdat, clvars ) 
    817854 
    818855      !----------------------------------------------------------------------- 
    … …  
    824861      DEALLOCATE( inpfiles ) 
    825862 
    826    END SUBROUTINE obs_rea_pro_dri 
     863   END SUBROUTINE obs_rea_prof 
    827864 
    828865END MODULE obs_read_prof 
Note: See TracChangeset for help on using the changeset viewer.