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 5659 for branches/2015 – NEMO

Changeset 5659 for branches/2015


Ignore:
Timestamp:
2015-07-31T11:59:15+02:00 (9 years ago)
Author:
mattmartin
Message:

Updated OBS simplification branch to the head of the trunk.

Location:
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
1 added
19 deleted
6 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r4990 r5659  
    2323   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch 
    2424   USE obs_read_prof            ! Reading and allocation of observations (Coriolis) 
    25    USE obs_read_sla             ! Reading and allocation of SLA observations   
    26    USE obs_read_sst             ! Reading and allocation of SST observations   
     25   USE obs_read_surf            ! Reading and allocation of SLA observations   
    2726   USE obs_readmdt              ! Reading and allocation of MDT for SLA. 
    28    USE obs_read_seaice          ! Reading and allocation of Sea Ice observations   
    29    USE obs_read_vel             ! Reading and allocation of velocity component observations 
    3027   USE obs_prep                 ! Preparation of obs. (grid search etc). 
    3128   USE obs_oper                 ! Observation operators 
     
    3431   USE obs_read_altbias         ! Bias treatment for altimeter 
    3532   USE obs_profiles_def         ! Profile data definitions 
    36    USE obs_profiles             ! Profile data storage 
    3733   USE obs_surf_def             ! Surface data definitions 
    38    USE obs_sla                  ! SLA data storage 
    39    USE obs_sst                  ! SST data storage 
    40    USE obs_seaice               ! Sea Ice data storage 
    4134   USE obs_types                ! Definitions for observation types 
    4235   USE mpp_map                  ! MPP mapping 
     
    6356   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles 
    6457   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles 
    65    LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set 
    66    LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set 
    67    LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles 
    6858   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies  
    69    LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files 
    70    LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files 
    7159   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature 
    72    LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature 
    73    LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data 
    74    LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files 
    7560   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration 
    7661   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations 
    77    LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data  
    78    LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data  
    79    LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data  
    80    LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data  
    81    LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files 
    8262   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height 
    8363   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity 
     
    9171   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS 
    9272   
     73   INTEGER, PUBLIC :: numobtypes   !: Number of observation types to read in. 
    9374   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method 
    9475   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method  
    95  
     76   INTEGER, DIMENSION(:), ALLOCATABLE :: nvarsprof !Number of profile variables 
     77   INTEGER, DIMENSION(:), ALLOCATABLE :: nextrprof !Number of profile extra variables 
     78   INTEGER, DIMENSION(:), ALLOCATABLE :: nvarssurf !Number of surface variables 
     79   INTEGER, DIMENSION(:), ALLOCATABLE :: nextrsurf !Number of surface extra variables 
    9680   INTEGER, DIMENSION(imaxavtypes) :: & 
    97       & endailyavtypes !: ENACT data types which are daily average 
    98  
     81      & dailyavtypes !: Data types which are daily average 
     82 
     83   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdata   ! Initial surface data 
     84   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdataqc ! Surface data after quality control 
     85   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdata   ! Initial profile data 
     86   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdataqc ! Profile data after quality control 
     87 
     88   CHARACTER(len=6),   PUBLIC, DIMENSION(:),   ALLOCATABLE :: obstypesprof 
     89   CHARACTER(len=6),   PUBLIC, DIMENSION(:),   ALLOCATABLE :: obstypessurf 
     90 
     91 
     92       
    9993   INTEGER, PARAMETER :: MaxNumFiles = 1000 
     94    
    10095   LOGICAL, DIMENSION(MaxNumFiles) :: & 
    10196      & ln_profb_ena, & !: Is the feedback files from ENACT data ? 
    102    !                    !: If so use endailyavtypes 
     97   !                    !: If so use dailyavtypes 
    10398      & ln_profb_enatim !: Change tim for 820 enact data set. 
    10499    
     
    106101      & ln_velfb_av   !: Is the velocity feedback files daily average? 
    107102   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    108       & ld_enact     !: Profile data is ENACT so use endailyavtypes 
     103      & ld_enact     !: Profile data is ENACT so use dailyavtypes 
    109104   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    110105      & ld_velav     !: Velocity data is daily averaged 
     
    135130      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    136131      !!        !  07-03  (K. Mogensen) General handling of profiles 
     132      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    137133      !!---------------------------------------------------------------------- 
    138134 
     
    140136 
    141137      !! * Local declarations 
    142       CHARACTER(len=128) :: enactfiles(MaxNumFiles) 
    143       CHARACTER(len=128) :: coriofiles(MaxNumFiles) 
    144138      CHARACTER(len=128) :: profbfiles(MaxNumFiles) 
    145       CHARACTER(len=128) :: sstfiles(MaxNumFiles)       
    146       CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)  
    147       CHARACTER(len=128) :: slafilesact(MaxNumFiles)       
    148       CHARACTER(len=128) :: slafilespas(MaxNumFiles)       
     139      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 
    149140      CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 
    150       CHARACTER(len=128) :: seaicefiles(MaxNumFiles)            
    151       CHARACTER(len=128) :: velcurfiles(MaxNumFiles)   
    152       CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)     
    153       CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 
    154       CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 
    155       CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 
    156       CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 
     141      CHARACTER(len=128) :: seaicefbfiles(MaxNumFiles) 
    157142      CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 
    158       CHARACTER(LEN=128) :: reysstname 
    159       CHARACTER(LEN=12)  :: reysstfmt 
    160143      CHARACTER(LEN=128) :: bias_file 
    161144      CHARACTER(LEN=20)  :: datestr=" ", timestr=" " 
    162       NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       & 
    163          &            ln_sla, ln_sladt, ln_slafb,                     & 
    164          &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       & 
    165          &            enactfiles, coriofiles, profbfiles,             & 
    166          &            slafilesact, slafilespas, slafbfiles,           & 
    167          &            sstfiles, sstfbfiles,                           & 
    168          &            ln_seaice, seaicefiles,                         & 
     145 
     146      NAMELIST/namobs/ln_t3d, ln_s3d, ln_sla, ln_sss, ln_ssh,         & 
     147         &            ln_sst, ln_seaice, ln_vel3d,                    & 
     148         &            ln_altbias, ln_nea, ln_grid_global,             & 
     149         &            ln_grid_search_lookup, ln_cl4,                  & 
     150         &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     151         &            ln_profb_ena, ln_profb_enatim,                  & 
     152         &            profbfiles, slafbfiles, sssfbfiles,             & 
     153         &            sshfbfiles, sstfbfiles, seaicefbfiles,          & 
     154         &            velfbfiles, bias_file, grid_search_file,        & 
    169155         &            dobsini, dobsend, n1dint, n2dint,               & 
    170156         &            nmsshc, mdtcorr, mdtcutoff,                     & 
    171          &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    & 
    172          &            ln_sstnight,                                    & 
    173          &            ln_grid_search_lookup,                          & 
    174          &            grid_search_file, grid_search_res,              & 
    175          &            ln_grid_global, bias_file, ln_altbias,          & 
    176          &            endailyavtypes, ln_s_at_t, ln_profb_ena,        & 
    177          &            ln_vel3d, ln_velavcur, velavcurfiles,           & 
    178          &            ln_velhrcur, velhrcurfiles,                     & 
    179          &            ln_velavadcp, velavadcpfiles,                   & 
    180          &            ln_velhradcp, velhradcpfiles,                   & 
    181          &            ln_velfb, velfbfiles, ln_velfb_av,              & 
    182          &            ln_profb_enatim, ln_ignmis, ln_cl4 
    183  
    184       INTEGER :: jprofset 
    185       INTEGER :: jveloset 
    186       INTEGER :: jvar 
    187       INTEGER :: jnumenact 
    188       INTEGER :: jnumcorio 
    189       INTEGER :: jnumprofb 
    190       INTEGER :: jnumslaact 
    191       INTEGER :: jnumslapas 
    192       INTEGER :: jnumslafb 
    193       INTEGER :: jnumsst 
    194       INTEGER :: jnumsstfb 
    195       INTEGER :: jnumseaice 
    196       INTEGER :: jnumvelavcur 
    197       INTEGER :: jnumvelhrcur   
    198       INTEGER :: jnumvelavadcp 
    199       INTEGER :: jnumvelhradcp    
    200       INTEGER :: jnumvelfb 
    201       INTEGER :: ji 
    202       INTEGER :: jset 
     157         &            grid_search_res, dailyavtypes 
     158 
     159      INTEGER :: jtype 
    203160      INTEGER :: ios                 ! Local integer output status for namelist read 
    204       LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 
    205  
     161      INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilesprof 
     162      INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilessurf 
     163      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilesprof 
     164      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilessurf 
     165      LOGICAL :: lmask(MaxNumFiles) 
    206166      !----------------------------------------------------------------------- 
    207167      ! Read namelist parameters 
    208168      !----------------------------------------------------------------------- 
    209169 
    210       enactfiles(:) = '' 
    211       coriofiles(:) = '' 
    212170      profbfiles(:) = '' 
    213       slafilesact(:) = '' 
    214       slafilespas(:) = '' 
    215171      slafbfiles(:) = '' 
    216       sstfiles(:)   = '' 
    217172      sstfbfiles(:) = '' 
    218       seaicefiles(:) = '' 
    219       velcurfiles(:) = '' 
    220       veladcpfiles(:) = '' 
    221       velavcurfiles(:) = '' 
    222       velhrcurfiles(:) = '' 
    223       velavadcpfiles(:) = '' 
    224       velhradcpfiles(:) = '' 
     173      seaicefbfiles(:) = '' 
    225174      velfbfiles(:) = '' 
    226       velcurfiles(:) = '' 
    227       veladcpfiles(:) = '' 
    228       endailyavtypes(:) = -1 
    229       endailyavtypes(1) = 820 
     175      dailyavtypes(:) = -1 
     176      dailyavtypes(1) = 820 
    230177      ln_profb_ena(:) = .FALSE. 
    231178      ln_profb_enatim(:) = .TRUE. 
    232179      ln_velfb_av(:) = .FALSE. 
    233180      ln_ignmis = .FALSE. 
    234        
     181 
    235182      CALL ini_date( dobsini ) 
    236183      CALL fin_date( dobsend ) 
    237   
     184 
    238185      ! Read Namelist namobs : control observation diagnostics 
    239186      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation 
     
    246193      IF(lwm) WRITE ( numond, namobs ) 
    247194 
    248       ! Count number of files for each type 
    249       IF (ln_ena) THEN 
    250          lmask(:) = .FALSE. 
    251          WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 
    252          jnumenact = COUNT(lmask) 
    253       ENDIF 
    254       IF (ln_cor) THEN 
    255          lmask(:) = .FALSE. 
    256          WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 
    257          jnumcorio = COUNT(lmask) 
    258       ENDIF 
    259       IF (ln_profb) THEN 
    260          lmask(:) = .FALSE. 
    261          WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 
    262          jnumprofb = COUNT(lmask) 
    263       ENDIF 
    264       IF (ln_sladt) THEN 
    265          lmask(:) = .FALSE. 
    266          WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 
    267          jnumslaact = COUNT(lmask) 
    268          lmask(:) = .FALSE. 
    269          WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 
    270          jnumslapas = COUNT(lmask) 
    271       ENDIF 
    272       IF (ln_slafb) THEN 
    273          lmask(:) = .FALSE. 
    274          WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 
    275          jnumslafb = COUNT(lmask) 
    276          lmask(:) = .FALSE. 
    277       ENDIF 
    278       IF (ln_ghrsst) THEN 
    279          lmask(:) = .FALSE. 
    280          WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 
    281          jnumsst = COUNT(lmask) 
    282       ENDIF       
    283       IF (ln_sstfb) THEN 
    284          lmask(:) = .FALSE. 
    285          WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 
    286          jnumsstfb = COUNT(lmask) 
    287          lmask(:) = .FALSE. 
    288       ENDIF 
    289       IF (ln_seaice) THEN 
    290          lmask(:) = .FALSE. 
    291          WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 
    292          jnumseaice = COUNT(lmask) 
    293       ENDIF 
    294       IF (ln_velavcur) THEN 
    295          lmask(:) = .FALSE. 
    296          WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 
    297          jnumvelavcur = COUNT(lmask) 
    298       ENDIF 
    299       IF (ln_velhrcur) THEN 
    300          lmask(:) = .FALSE. 
    301          WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 
    302          jnumvelhrcur = COUNT(lmask) 
    303       ENDIF 
    304       IF (ln_velavadcp) THEN 
    305          lmask(:) = .FALSE. 
    306          WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 
    307          jnumvelavadcp = COUNT(lmask) 
    308       ENDIF 
    309       IF (ln_velhradcp) THEN 
    310          lmask(:) = .FALSE. 
    311          WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 
    312          jnumvelhradcp = COUNT(lmask) 
    313       ENDIF 
    314       IF (ln_velfb) THEN 
    315          lmask(:) = .FALSE. 
    316          WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 
    317          jnumvelfb = COUNT(lmask) 
    318          lmask(:) = .FALSE. 
    319       ENDIF 
    320        
    321       ! Control print 
     195      !Set up list of observation types to be used 
     196      numproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     197      numsurftypes = COUNT( (/ln_sla, ln_sss, ln_sst, ln_seaice /) ) 
     198      IF ( numproftypes > 0 ) THEN 
     199       
     200         ALLOCATE( obstypesprof(numproftypes) ) 
     201         ALLOCATE( jnumfilesprof(numproftypes) ) 
     202         ALLOCATE( obsfilesprof(numproftypes, MaxNumFiles) ) 
     203       
     204         DO jtype = 1, numproftypes 
     205            IF (ln_t3d .OR. ln_s3d) THEN 
     206               obsfilesprof(:,jtype) = profbfiles(:) 
     207               obstypesprof(jtype) = 'prof' 
     208            ENDIF 
     209            IF (ln_vel3d) THEN 
     210               obsfilesprof(:,jtype) = velfbfiles(:) 
     211               obstypesprof(jtype) = 'vel' 
     212            ENDIF 
     213          
     214            lmask(:) = .FALSE. 
     215            WHERE (obsfilesprof(jtype,:) /= '') lmask(:) = .TRUE. 
     216            jnumfilesprof(jtype) = COUNT(lmask) 
     217         END DO 
     218          
     219      ENDIF 
     220       
     221      IF ( numsurftypes > 0 ) THEN 
     222       
     223         ALLOCATE( obstypessurf(numsurftypes) ) 
     224         ALLOCATE( jnumfilessurf(numproftypes) ) 
     225         ALLOCATE( obsfilessurf(numsurftypes, MaxNumFiles) ) 
     226          
     227         DO jtype = 1, numsurftypes 
     228            IF (ln_sla) THEN 
     229               obsfilessurf(:,jtype) = slafbfiles(:) 
     230               obstypessurf(jtype) = 'sla' 
     231            ENDIF 
     232            IF (ln_sss) THEN 
     233               obsfilessurf(:,jtype) = sssfbfiles(:) 
     234               obstypessurf(jtype) = 'sss' 
     235            ENDIF 
     236            IF (ln_sst) THEN 
     237               obsfilessurf(:,jtype) = sstfbfiles(:) 
     238               obstypessurf(jtype) = 'sst' 
     239            ENDIF 
     240#if defined key_lim2 || defined key_lim3 
     241            IF (ln_seaice) THEN 
     242               obsfilessurf(:,jtype) = seaicefbfiles(:) 
     243               obstypessurf(jtype) = 'seaice' 
     244            ENDIF 
     245#endif 
     246 
     247            lmask(:) = .FALSE. 
     248            WHERE (obsfilessurf(jtype,:) /= '') lmask(:) = .TRUE. 
     249            jnumfilessurf(jtype) = COUNT(lmask) 
     250          
     251         END DO 
     252          
     253      ENDIF 
     254 
     255      !Write namelist settings to stdout 
    322256      IF(lwp) THEN 
    323257         WRITE(numout,*) 
     
    327261         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d 
    328262         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d 
    329          WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena 
    330          WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor 
    331          WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb 
    332263         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla 
    333          WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt 
    334          WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb 
    335264         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh 
    336265         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst 
    337          WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst     
    338          WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst 
    339          WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb 
    340266         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight 
    341267         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss 
    342268         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice 
    343269         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d 
    344          WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur 
    345          WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur 
    346          WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp 
    347          WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp 
    348          WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb 
    349          WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global 
     270         WRITE(numout,*) '             Global distribution of observations        ln_grid_global = ',ln_grid_global 
    350271         WRITE(numout,*) & 
    351272   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup 
    352273         IF (ln_grid_search_lookup) & 
    353274            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file 
    354          IF (ln_ena) THEN 
    355             DO ji = 1, jnumenact 
    356                WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', & 
    357                   TRIM(enactfiles(ji)) 
    358             END DO 
    359          ENDIF 
    360          IF (ln_cor) THEN 
    361             DO ji = 1, jnumcorio 
    362                WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', & 
    363                   TRIM(coriofiles(ji)) 
    364             END DO 
    365          ENDIF 
    366          IF (ln_profb) THEN 
    367             DO ji = 1, jnumprofb 
    368                IF (ln_profb_ena(ji)) THEN 
    369                   WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', & 
    370                      TRIM(profbfiles(ji)) 
    371                ELSE 
    372                   WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', & 
    373                      TRIM(profbfiles(ji)) 
    374                ENDIF 
    375                WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji) 
    376             END DO 
    377          ENDIF 
    378          IF (ln_sladt) THEN 
    379             DO ji = 1, jnumslaact 
    380                WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', & 
    381                   TRIM(slafilesact(ji)) 
    382             END DO 
    383             DO ji = 1, jnumslapas 
    384                WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', & 
    385                   TRIM(slafilespas(ji)) 
    386             END DO 
    387          ENDIF 
    388          IF (ln_slafb) THEN 
    389             DO ji = 1, jnumslafb 
    390                WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', & 
    391                   TRIM(slafbfiles(ji)) 
    392             END DO 
    393          ENDIF 
    394          IF (ln_ghrsst) THEN 
    395             DO ji = 1, jnumsst 
    396                WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', & 
    397                   TRIM(sstfiles(ji)) 
    398             END DO 
    399          ENDIF 
    400          IF (ln_sstfb) THEN 
    401             DO ji = 1, jnumsstfb 
    402                WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', & 
    403                   TRIM(sstfbfiles(ji)) 
    404             END DO 
    405          ENDIF 
    406          IF (ln_seaice) THEN 
    407             DO ji = 1, jnumseaice 
    408                WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', & 
    409                   TRIM(seaicefiles(ji)) 
    410             END DO 
    411          ENDIF 
    412          IF (ln_velavcur) THEN 
    413             DO ji = 1, jnumvelavcur 
    414                WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', & 
    415                   TRIM(velavcurfiles(ji)) 
    416             END DO 
    417          ENDIF 
    418          IF (ln_velhrcur) THEN 
    419             DO ji = 1, jnumvelhrcur 
    420                WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', & 
    421                   TRIM(velhrcurfiles(ji)) 
    422             END DO 
    423          ENDIF 
    424          IF (ln_velavadcp) THEN 
    425             DO ji = 1, jnumvelavadcp 
    426                WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', & 
    427                   TRIM(velavadcpfiles(ji)) 
    428             END DO 
    429          ENDIF 
    430          IF (ln_velhradcp) THEN 
    431             DO ji = 1, jnumvelhradcp 
    432                WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', & 
    433                   TRIM(velhradcpfiles(ji)) 
    434             END DO 
    435          ENDIF 
    436          IF (ln_velfb) THEN 
    437             DO ji = 1, jnumvelfb 
    438                IF (ln_velfb_av(ji)) THEN 
    439                   WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', & 
    440                      TRIM(velfbfiles(ji)) 
    441                ELSE 
    442                   WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', & 
    443                      TRIM(velfbfiles(ji)) 
    444                ENDIF 
    445             END DO 
    446          ENDIF 
    447          WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini 
     275         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsin 
    448276         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend 
    449277         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint 
     
    455283         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias 
    456284         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis 
    457          WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes 
     285         WRITE(numout,*) '             Daily average types                                   = ', dailyavtypes 
     286 
     287         IF ( numproftypes > 0 ) THEN 
     288            DO jtype = 1, numproftypes 
     289               DO ji = 1, jnumfilesprof(jtype) 
     290                  WRITE(numout,'(1X,2A)') '             '//obstypesprof(jtype)//' input observation file names  = ', & 
     291                     TRIM(obsfilesprof(jtype,ji)) 
     292                  IF ( TRIM(obstypesprof(jtype)) == 'prof' ) & 
     293                     WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji) 
     294               END DO 
     295            END DO 
     296         ENDIF 
     297          
     298         IF ( numsurftypes > 0 ) THEN 
     299            DO jtype = 1, numsurftypes 
     300               DO ji = 1, jnumfilessurf(jtype) 
     301                  WRITE(numout,'(1X,2A)') '             '//obstypessurf(jtype)//' input observation file names  = ', & 
     302                     TRIM(obsfilessurf(jtype,ji)) 
     303               END DO 
     304            END DO 
     305         ENDIF 
    458306 
    459307      ENDIF 
     
    470318      ! Parameter control 
    471319#if defined key_diaobs 
    472       IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 
    473          & ( .NOT. ln_vel3d ).AND.                                         & 
    474          & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 
    475          & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 
     320      IF ( numobtypes == 0 ) THEN 
    476321         IF(lwp) WRITE(numout,cform_war) 
    477322         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 
     
    494339      ! Depending on switches read the various observation types 
    495340      !----------------------------------------------------------------------- 
    496       !  - Temperature/salinity profiles 
    497  
    498       IF ( ln_t3d .OR. ln_s3d ) THEN 
    499  
    500          ! Set the number of variables for profiles to 2 (T and S) 
    501          nprofvars = 2 
    502          ! Set the number of extra variables for profiles to 1 (insitu temp). 
    503          nprofextr = 1 
    504  
    505          ! Count how may insitu data sets we have and allocate data. 
    506          jprofset = 0 
    507          IF ( ln_ena ) jprofset = jprofset + 1 
    508          IF ( ln_cor ) jprofset = jprofset + 1 
    509          IF ( ln_profb ) jprofset = jprofset + jnumprofb 
    510          nprofsets = jprofset 
    511          IF ( nprofsets > 0 ) THEN 
    512             ALLOCATE(ld_enact(nprofsets)) 
    513             ALLOCATE(profdata(nprofsets)) 
    514             ALLOCATE(prodatqc(nprofsets)) 
    515          ENDIF 
    516  
    517          jprofset = 0 
    518            
    519          ! ENACT insitu data 
    520  
    521          IF ( ln_ena ) THEN 
    522  
    523             jprofset = jprofset + 1 
     341       
     342      IF ( numproftypes > 0 ) THEN 
     343       
     344         ALLOCATE(profdata(numproftypes)) 
     345         ALLOCATE(profdataqc(numproftypes)) 
     346         ALLOCATE(nvarsprof(numproftypes)) 
     347         ALLOCATE(nextrprof(numproftypes)) 
    524348             
    525             ld_enact(jprofset) = .TRUE. 
    526  
    527             CALL obs_rea_pro_dri( 1, profdata(jprofset),          & 
    528                &                  jnumenact, enactfiles(1:jnumenact), & 
    529                &                  nprofvars, nprofextr,        & 
    530                &                  nitend-nit000+2,             & 
    531                &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    532                &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 
    533                &                  kdailyavtypes = endailyavtypes ) 
    534  
    535             DO jvar = 1, 2 
    536  
    537                CALL obs_prof_staend( profdata(jprofset), jvar ) 
    538  
     349         DO jtype = 1, numproftypes 
     350       
     351            nvarsprof(jtype) = 2 
     352            IF ( TRIM(obstypesprof(jtype)) == 'prof' ) nextrprof(jtype) = 1 
     353            IF ( TRIM(obstypesprof(jtype)) == 'vel' )  nextrprof(jtype) = 2 
     354 
     355            !Read in profile or velocity obs types 
     356            CALL obs_rea_prof( profdata(jtype),          & 
     357               &               jnumfilesprof(jtype),       & 
     358               &               obsfilesprof(jtype,1:jnumfilesprof(jtype)), & 
     359               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2,             & 
     360               &               dobsini, dobsend, ln_t3d, ln_s3d, & 
     361               &               ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 
     362               &               kdailyavtypes = dailyavtypes ) 
     363             
     364            DO jvar = 1, nvars 
     365               CALL obs_prof_staend( profdata(jtype), jvar ) 
    539366            END DO 
    540  
    541             CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
     367          
     368            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype),   & 
    542369               &              ln_t3d, ln_s3d, ln_nea, & 
    543                &              kdailyavtypes=endailyavtypes ) 
    544              
    545          ENDIF 
    546  
    547          ! Coriolis insitu data 
    548  
    549          IF ( ln_cor ) THEN 
    550             
    551             jprofset = jprofset + 1 
    552  
    553             ld_enact(jprofset) = .FALSE. 
    554  
    555             CALL obs_rea_pro_dri( 2, profdata(jprofset),          & 
    556                &                  jnumcorio, coriofiles(1:jnumcorio), & 
    557                &                  nprofvars, nprofextr,        & 
    558                &                  nitend-nit000+2,             & 
    559                &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    560                &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 
    561  
    562             DO jvar = 1, 2 
    563  
    564                CALL obs_prof_staend( profdata(jprofset), jvar ) 
    565  
    566             END DO 
    567  
    568             CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    569                  &            ln_t3d, ln_s3d, ln_nea ) 
    570              
    571          ENDIF 
    572   
    573          ! Feedback insitu data 
    574  
    575          IF ( ln_profb ) THEN 
    576             
    577             DO jset = 1, jnumprofb 
    578                 
    579                jprofset = jprofset + 1 
    580                ld_enact (jprofset) = ln_profb_ena(jset) 
    581  
    582                CALL obs_rea_pro_dri( 0, profdata(jprofset),          & 
    583                   &                  1, profbfiles(jset:jset), & 
    584                   &                  nprofvars, nprofextr,        & 
    585                   &                  nitend-nit000+2,             & 
    586                   &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    587                   &                  ln_ignmis, ln_s_at_t, & 
    588                   &                  ld_enact(jprofset).AND.& 
    589                   &                  ln_profb_enatim(jset), & 
    590                   &                  .FALSE., kdailyavtypes = endailyavtypes ) 
    591                 
    592                DO jvar = 1, 2 
    593                    
    594                   CALL obs_prof_staend( profdata(jprofset), jvar ) 
    595                    
    596                END DO 
    597                 
    598                IF ( ld_enact(jprofset) ) THEN 
    599                   CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    600                      &              ln_t3d, ln_s3d, ln_nea, & 
    601                      &              kdailyavtypes = endailyavtypes ) 
    602                ELSE 
    603                   CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    604                      &              ln_t3d, ln_s3d, ln_nea ) 
    605                ENDIF 
    606                 
    607             END DO 
    608  
    609          ENDIF 
    610  
    611       ENDIF 
    612  
    613       !  - Sea level anomalies 
    614       IF ( ln_sla ) THEN 
    615         ! Set the number of variables for sla to 1 
    616          nslavars = 1 
    617  
    618          ! Set the number of extra variables for sla to 2 
    619          nslaextr = 2 
    620           
    621          ! Set the number of sla data sets to 2 
    622          nslasets = 0 
    623          IF ( ln_sladt ) THEN 
    624             nslasets = nslasets + 2 
    625          ENDIF 
    626          IF ( ln_slafb ) THEN 
    627             nslasets = nslasets + jnumslafb 
    628          ENDIF 
    629           
    630          ALLOCATE(sladata(nslasets)) 
    631          ALLOCATE(sladatqc(nslasets)) 
    632          sladata(:)%nsurf=0 
    633          sladatqc(:)%nsurf=0 
    634  
    635          nslasets = 0 
    636  
    637          ! AVISO SLA data 
    638  
    639          IF ( ln_sladt ) THEN 
    640  
    641             ! Active SLA observations 
    642              
    643             nslasets = nslasets + 1 
    644              
    645             CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 
    646                &              slafilesact(1:jnumslaact), & 
    647                &              nslavars, nslaextr, nitend-nit000+2, & 
    648                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    649             CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    650                &              ln_sla, ln_nea ) 
    651              
    652             ! Passive SLA observations 
    653              
    654             nslasets = nslasets + 1 
    655              
    656             CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 
    657                &              slafilespas(1:jnumslapas), & 
    658                &              nslavars, nslaextr, nitend-nit000+2, & 
    659                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    660              
    661             CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    662                &              ln_sla, ln_nea ) 
    663  
    664          ENDIF 
    665           
    666          ! Feedback SLA data 
    667  
    668          IF ( ln_slafb ) THEN 
    669  
    670             DO jset = 1, jnumslafb 
    671              
    672                nslasets = nslasets + 1 
    673              
    674                CALL obs_rea_sla( 0, sladata(nslasets), 1, & 
    675                   &              slafbfiles(jset:jset), & 
    676                   &              nslavars, nslaextr, nitend-nit000+2, & 
    677                   &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    678                CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    679                   &              ln_sla, ln_nea ) 
    680  
    681             END DO                
    682  
    683          ENDIF 
    684           
    685          CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 
    686              
    687          ! read in altimeter bias 
    688           
    689          IF ( ln_altbias ) THEN      
    690             CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 
    691          ENDIF 
    692       
    693       ENDIF 
    694  
    695       !  - Sea surface height 
    696       IF ( ln_ssh ) THEN 
    697          IF(lwp) WRITE(numout,*) ' SSH currently not available' 
    698       ENDIF 
    699  
    700       !  - Sea surface temperature 
    701       IF ( ln_sst ) THEN 
    702  
    703          ! Set the number of variables for sst to 1 
    704          nsstvars = 1 
    705  
    706          ! Set the number of extra variables for sst to 0 
    707          nsstextr = 0 
    708  
    709          nsstsets = 0 
    710  
    711          IF (ln_reysst) nsstsets = nsstsets + 1 
    712          IF (ln_ghrsst) nsstsets = nsstsets + 1 
    713          IF ( ln_sstfb ) THEN 
    714             nsstsets = nsstsets + jnumsstfb 
    715          ENDIF 
    716  
    717          ALLOCATE(sstdata(nsstsets)) 
    718          ALLOCATE(sstdatqc(nsstsets)) 
    719          ALLOCATE(ld_sstnight(nsstsets)) 
    720          sstdata(:)%nsurf=0 
    721          sstdatqc(:)%nsurf=0     
    722          ld_sstnight(:)=.false. 
    723  
    724          nsstsets = 0 
    725  
    726          IF (ln_reysst) THEN 
    727  
    728             nsstsets = nsstsets + 1 
    729  
    730             ld_sstnight(nsstsets) = ln_sstnight 
    731  
    732             CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 
    733                &                  nsstvars, nsstextr, & 
    734                &                  nitend-nit000+2, dobsini, dobsend ) 
    735             CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 
    736                &              ln_nea ) 
    737  
    738         ENDIF 
    739          
    740         IF (ln_ghrsst) THEN 
    741          
    742             nsstsets = nsstsets + 1 
    743  
    744             ld_sstnight(nsstsets) = ln_sstnight 
    745            
    746             CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 
    747                &              sstfiles(1:jnumsst), & 
    748                &              nsstvars, nsstextr, nitend-nit000+2, & 
    749                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    750             CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 
    751                &              ln_nea ) 
    752  
    753         ENDIF 
    754                 
    755          ! Feedback SST data 
    756  
    757          IF ( ln_sstfb ) THEN 
    758  
    759             DO jset = 1, jnumsstfb 
    760              
    761                nsstsets = nsstsets + 1 
    762  
    763                ld_sstnight(nsstsets) = ln_sstnight 
    764              
    765                CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 
    766                   &              sstfbfiles(jset:jset), & 
    767                   &              nsstvars, nsstextr, nitend-nit000+2, & 
    768                   &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    769                CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 
    770                   &              ln_sst, ln_nea ) 
    771  
    772             END DO                
    773  
    774          ENDIF 
    775  
    776       ENDIF 
    777  
    778       !  - Sea surface salinity 
    779       IF ( ln_sss ) THEN 
    780          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    781       ENDIF 
    782  
    783       !  - Sea Ice Concentration 
    784        
    785       IF ( ln_seaice ) THEN 
    786  
    787          ! Set the number of variables for seaice to 1 
    788          nseaicevars = 1 
    789  
    790          ! Set the number of extra variables for seaice to 0 
    791          nseaiceextr = 0 
    792           
    793          ! Set the number of data sets to 1 
    794          nseaicesets = 1 
    795  
    796          ALLOCATE(seaicedata(nseaicesets)) 
    797          ALLOCATE(seaicedatqc(nseaicesets)) 
    798          seaicedata(:)%nsurf=0 
    799          seaicedatqc(:)%nsurf=0 
    800  
    801          CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 
    802             &                 seaicefiles(1:jnumseaice), & 
    803             &                 nseaicevars, nseaiceextr, nitend-nit000+2, & 
    804             &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
    805  
    806          CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 
    807             &                 ln_seaice, ln_nea ) 
    808   
    809       ENDIF 
    810  
    811       IF (ln_vel3d) THEN 
    812  
    813          ! Set the number of variables for profiles to 2 (U and V) 
    814          nvelovars = 2 
    815  
    816          ! Set the number of extra variables for profiles to 2 to store  
    817          ! rotation parameters 
    818          nveloextr = 2 
    819  
    820          jveloset = 0 
    821           
    822          IF ( ln_velavcur ) jveloset = jveloset + 1 
    823          IF ( ln_velhrcur ) jveloset = jveloset + 1 
    824          IF ( ln_velavadcp ) jveloset = jveloset + 1 
    825          IF ( ln_velhradcp ) jveloset = jveloset + 1 
    826          IF (ln_velfb) jveloset = jveloset + jnumvelfb 
    827  
    828          nvelosets = jveloset 
    829          IF ( nvelosets > 0 ) THEN 
    830             ALLOCATE( velodata(nvelosets) ) 
    831             ALLOCATE( veldatqc(nvelosets) ) 
    832             ALLOCATE( ld_velav(nvelosets) ) 
    833          ENDIF 
    834           
    835          jveloset = 0 
    836           
    837          ! Daily averaged data 
    838  
    839          IF ( ln_velavcur ) THEN 
    840              
    841             jveloset = jveloset + 1 
    842              
    843             ld_velav(jveloset) = .TRUE. 
    844              
    845             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 
    846                &                  velavcurfiles(1:jnumvelavcur), & 
    847                &                  nvelovars, nveloextr, & 
    848                &                  nitend-nit000+2,              & 
    849                &                  dobsini, dobsend, ln_ignmis, & 
    850                &                  ld_velav(jveloset), & 
    851                &                  .FALSE. ) 
    852              
    853             DO jvar = 1, 2 
    854                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    855             END DO 
    856              
    857             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    858                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    859              
    860          ENDIF 
    861  
    862          ! High frequency data 
    863  
    864          IF ( ln_velhrcur ) THEN 
    865              
    866             jveloset = jveloset + 1 
    867              
    868             ld_velav(jveloset) = .FALSE. 
    869                 
    870             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 
    871                &                  velhrcurfiles(1:jnumvelhrcur), & 
    872                &                  nvelovars, nveloextr, & 
    873                &                  nitend-nit000+2,              & 
    874                &                  dobsini, dobsend, ln_ignmis, & 
    875                &                  ld_velav(jveloset), & 
    876                &                  .FALSE. ) 
    877              
    878             DO jvar = 1, 2 
    879                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    880             END DO 
    881              
    882             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    883                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    884              
    885          ENDIF 
    886  
    887          ! Daily averaged data 
    888  
    889          IF ( ln_velavadcp ) THEN 
    890              
    891             jveloset = jveloset + 1 
    892              
    893             ld_velav(jveloset) = .TRUE. 
    894              
    895             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 
    896                &                  velavadcpfiles(1:jnumvelavadcp), & 
    897                &                  nvelovars, nveloextr, & 
    898                &                  nitend-nit000+2,              & 
    899                &                  dobsini, dobsend, ln_ignmis, & 
    900                &                  ld_velav(jveloset), & 
    901                &                  .FALSE. ) 
    902              
    903             DO jvar = 1, 2 
    904                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    905             END DO 
    906              
    907             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    908                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    909              
    910          ENDIF 
    911  
    912          ! High frequency data 
    913  
    914          IF ( ln_velhradcp ) THEN 
    915              
    916             jveloset = jveloset + 1 
    917              
    918             ld_velav(jveloset) = .FALSE. 
    919                 
    920             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 
    921                &                  velhradcpfiles(1:jnumvelhradcp), & 
    922                &                  nvelovars, nveloextr, & 
    923                &                  nitend-nit000+2,              & 
    924                &                  dobsini, dobsend, ln_ignmis, & 
    925                &                  ld_velav(jveloset), & 
    926                &                  .FALSE. ) 
    927              
    928             DO jvar = 1, 2 
    929                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    930             END DO 
    931              
    932             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    933                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    934              
    935          ENDIF 
    936  
    937          IF ( ln_velfb ) THEN 
    938  
    939             DO jset = 1, jnumvelfb 
    940              
    941                jveloset = jveloset + 1 
    942  
    943                ld_velav(jveloset) = ln_velfb_av(jset) 
    944                 
    945                CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 
    946                   &                  velfbfiles(jset:jset), & 
    947                   &                  nvelovars, nveloextr, & 
    948                   &                  nitend-nit000+2,              & 
    949                   &                  dobsini, dobsend, ln_ignmis, & 
    950                   &                  ld_velav(jveloset), & 
    951                   &                  .FALSE. ) 
    952                 
    953                DO jvar = 1, 2 
    954                   CALL obs_prof_staend( velodata(jveloset), jvar ) 
    955                END DO 
    956                 
    957                CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    958                   &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    959  
    960  
    961             END DO 
    962              
    963          ENDIF 
    964  
    965       ENDIF 
    966       
     370               &              kdailyavtypes = dailyavtypes ) 
     371 
     372         END DO 
     373 
     374         DEALLOCATE( jnumfilesprof, obsfilesprof ) 
     375 
     376      ENDIF 
     377 
     378      IF ( numsurftypes > 0 ) THEN 
     379       
     380         ALLOCATE(surfdata(numsurftypes)) 
     381         ALLOCATE(surfdatatqc(numsurftypes)) 
     382         ALLOCATE(nvarssurf(numsurftypes)) 
     383         ALLOCATE(nextrsurf(numsurftypes)) 
     384          
     385         DO jtype = 1, numsurftypes 
     386       
     387            nvarssurf(jtype) = 1 
     388            nextrsurf(jtype) = 0 
     389            IF ( TRIM(obstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
     390 
     391            !Read in surface obs types 
     392            CALL obs_rea_surf( surfdata(jtype), jnumfilessurf(jtype), & 
     393               &               obsfilessurf(jtype,1:jnumfilessurf(jtype)), & 
     394               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
     395               &               dobsini, dobsend, ln_ignmis, .FALSE. ) 
     396 
     397            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     398          
     399            IF ( TRIM(obstypessurf(jtype)) == 'sla' ) THEN 
     400               CALL obs_rea_mdt( surfdataqc(jtype), n2dint ) 
     401               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), n2dint, bias_file ) 
     402            ENDIF 
     403 
     404            DEALLOCATE( jnumfilessurf, obsfilessurf ) 
     405 
     406         END DO 
     407 
    967408   END SUBROUTINE dia_obs_init 
    968409 
     
    1017458      !! * Local declarations 
    1018459      INTEGER :: idaystp                ! Number of timesteps per day 
    1019       INTEGER :: jprofset               ! Profile data set loop variable 
    1020       INTEGER :: jslaset                ! SLA data set loop variable 
    1021       INTEGER :: jsstset                ! SST data set loop variable 
    1022       INTEGER :: jseaiceset             ! sea ice data set loop variable 
    1023       INTEGER :: jveloset               ! velocity profile data loop variable 
     460      INTEGER :: jtype                  ! data loop variable 
    1024461      INTEGER :: jvar                   ! Variable number     
    1025462#if ! defined key_lim2 && ! defined key_lim3 
     
    1050487      !----------------------------------------------------------------------- 
    1051488 
    1052       !  - Temperature/salinity profiles 
    1053       IF ( ln_t3d .OR. ln_s3d ) THEN 
    1054          DO jprofset = 1, nprofsets 
    1055             IF ( ld_enact(jprofset) ) THEN 
    1056                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1057                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1058                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1059                   &              gdept_1d, tmask, n1dint, n2dint,        & 
    1060                   &              kdailyavtypes = endailyavtypes ) 
    1061             ELSE 
    1062                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1063                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1064                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1065                   &              gdept_1d, tmask, n1dint, n2dint              ) 
    1066             ENDIF 
     489      IF ( numproftypes > 0 ) THEN 
     490         DO jtype = 1, numproftypes 
     491       
     492            SELECT CASE ( TRIM(obstypesprof(jtype)) ) 
     493            CASE('prof') 
     494               profvar1(:,:,:) = tsn(:,:,:,jp_tem) 
     495               profvar2(:,:,:) = tsn(:,:,:,jp_sal) 
     496               profmask1(:,:,:) = tmask(:,:,:) 
     497               profmask2(:,:,:) = tmask(:,:,:) 
     498            CASE('vel') 
     499               profvar1(:,:,:) = un(:,:,:) 
     500               profvar2(:,:,:) = vn(:,:,:) 
     501               profmask1(:,:,:) = umask(:,:,:) 
     502               profmask2(:,:,:) = vmask(:,:,:) 
     503            END SELECT 
     504             
     505            CALL obs_prof_opt( profdataqc(jtype),                       & 
     506               &               kstp, jpi, jpj, jpk, nit000, idaystp,   & 
     507               &               profvar1, profvar2,   & 
     508               &               gdept_1d, profmask1, profmask2, n1dint, n2dint,        & 
     509               &               kdailyavtypes = dailyavtypes ) 
     510             
    1067511         END DO 
    1068       ENDIF 
    1069  
    1070       !  - Sea surface anomaly 
    1071       IF ( ln_sla ) THEN 
    1072          DO jslaset = 1, nslasets 
    1073             CALL obs_sla_opt( sladatqc(jslaset),            & 
    1074                &              kstp, jpi, jpj, nit000, sshn, & 
    1075                &              tmask(:,:,1), n2dint ) 
    1076          END DO          
    1077       ENDIF 
    1078  
    1079       !  - Sea surface temperature 
    1080       IF ( ln_sst ) THEN 
    1081          DO jsstset = 1, nsstsets 
    1082             CALL obs_sst_opt( sstdatqc(jsstset),                & 
    1083                &              kstp, jpi, jpj, nit000, idaystp,  & 
    1084                &              tsn(:,:,1,jp_tem), tmask(:,:,1),  & 
    1085                &              n2dint, ld_sstnight(jsstset) ) 
     512 
     513      ENDIF 
     514 
     515      IF ( numsurftypes > 0 ) THEN 
     516         DO jtype = 1, numsurftypes 
     517       
     518            SELECT CASE ( TRIM(obstypessurf(jtype)) ) 
     519            CASE('sst') 
     520               surfvar(:,:) = tsn(:,:,1,jp_tem) 
     521            CASE('sla') 
     522               surfvar(:,:) = sshn(:,:) 
     523            CASE('sss') 
     524               surfvar(:,:) = tsn(:,:,1,jp_sal) 
     525#if defined key_lim2 || defined key_lim3 
     526            CASE('seaice') 
     527               surfvar(:,:) = 1._wp - frld(:,:) 
     528#endif 
     529            END SELECT 
     530          
     531            CALL obs_surf_opt( surfdatqc(jtype),             & 
     532               &               kstp, jpi, jpj, nit000, surfvar, & 
     533               &               tmask(:,:,1), n2dint, ld_sstnight ) 
     534             
    1086535         END DO 
    1087       ENDIF 
    1088  
    1089       !  - Sea surface salinity 
    1090       IF ( ln_sss ) THEN 
    1091          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    1092       ENDIF 
    1093  
    1094 #if defined key_lim2 || defined key_lim3 
    1095       IF ( ln_seaice ) THEN 
    1096          DO jseaiceset = 1, nseaicesets 
    1097             CALL obs_seaice_opt( seaicedatqc(jseaiceset),      & 
    1098                &              kstp, jpi, jpj, nit000, 1.-frld, & 
    1099                &              tmask(:,:,1), n2dint ) 
    1100          END DO 
    1101       ENDIF       
    1102 #endif 
    1103  
    1104       !  - Velocity profiles 
    1105       IF ( ln_vel3d ) THEN 
    1106          DO jveloset = 1, nvelosets 
    1107            ! zonal component of velocity 
    1108            CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, & 
    1109               &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, & 
    1110                              n1dint, n2dint, ld_velav(jveloset) ) 
    1111          END DO 
    1112       ENDIF 
    1113  
     536          
     537      ENDIF 
     538       
    1114539#if ! defined key_lim2 && ! defined key_lim3 
    1115540      CALL wrk_dealloc(jpi,jpj,frld)  
     
    1139564      !! * Local declarations 
    1140565 
    1141       INTEGER :: jprofset                 ! Profile data set loop variable 
    1142       INTEGER :: jveloset                 ! Velocity data set loop variable 
    1143       INTEGER :: jslaset                  ! SLA data set loop variable 
    1144       INTEGER :: jsstset                  ! SST data set loop variable 
    1145       INTEGER :: jseaiceset               ! Sea Ice data set loop variable 
    1146       INTEGER :: jset 
    1147       INTEGER :: jfbini 
    1148       CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    1149       CHARACTER(LEN=10) :: cdtmp 
     566      INTEGER :: jtype                    ! Data set loop variable 
    1150567      !----------------------------------------------------------------------- 
    1151568      ! Depending on switches call various observation output routines 
    1152569      !----------------------------------------------------------------------- 
    1153570 
    1154       !  - Temperature/salinity profiles 
    1155  
    1156       IF( ln_t3d .OR. ln_s3d ) THEN 
    1157  
    1158          ! Copy data from prodatqc to profdata structures 
    1159          DO jprofset = 1, nprofsets 
    1160  
    1161             CALL obs_prof_decompress( prodatqc(jprofset), & 
    1162                  &                    profdata(jprofset), .TRUE., numout ) 
    1163  
     571      IF ( numproftypes > 0 ) THEN 
     572         DO jtype = 1, numproftypes 
     573       
     574            CALL obs_prof_decompress( profdataqc(jtype), & 
     575               &                      profdata(jtype), .TRUE., numout ) 
     576 
     577            CALL obs_wri_prof( obstypesprof(jtype), profdata(jtype), n2dint ) 
     578          
    1164579         END DO 
    1165  
    1166          ! Write the profiles. 
    1167  
    1168          jprofset = 0 
    1169  
    1170          ! ENACT insitu data 
    1171  
    1172          IF ( ln_ena ) THEN 
    1173             
    1174             jprofset = jprofset + 1 
    1175  
    1176             CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 
    1177  
    1178          ENDIF 
    1179  
    1180          ! Coriolis insitu data 
    1181  
    1182          IF ( ln_cor ) THEN 
    1183              
    1184             jprofset = jprofset + 1 
    1185  
    1186             CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 
    1187              
    1188          ENDIF 
    1189           
    1190          ! Feedback insitu data 
    1191  
    1192          IF ( ln_profb ) THEN 
    1193  
    1194             jfbini = jprofset + 1 
    1195  
    1196             DO jprofset = jfbini, nprofsets 
    1197                 
    1198                jset = jprofset - jfbini + 1 
    1199                WRITE(cdtmp,'(A,I2.2)')'profb_',jset 
    1200                CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 
    1201  
    1202             END DO 
    1203  
    1204          ENDIF 
    1205  
    1206       ENDIF 
    1207  
    1208       !  - Sea surface anomaly 
    1209       IF ( ln_sla ) THEN 
    1210  
    1211          ! Copy data from sladatqc to sladata structures 
    1212          DO jslaset = 1, nslasets 
    1213  
    1214               CALL obs_surf_decompress( sladatqc(jslaset), & 
    1215                  &                    sladata(jslaset), .TRUE., numout ) 
     580          
     581      ENDIF 
     582 
     583      IF ( numsurftypes > 0 ) THEN 
     584         DO jtype = 1, numsurftypes 
     585       
     586            CALL obs_surf_decompress( surfdatqc(jtype), & 
     587               &                      surfdata(jtype), .TRUE., numout ) 
     588 
     589            CALL obs_wri_surf( obstypessurf(jtype), surfdata(jtype), n2dint ) 
    1216590 
    1217591         END DO 
    1218  
    1219          jslaset = 0  
    1220  
    1221          ! Write the AVISO SLA data 
    1222  
    1223          IF ( ln_sladt ) THEN 
    1224              
    1225             jslaset = 1 
    1226             CALL obs_wri_sla( 'aviso_act', sladata(jslaset) ) 
    1227             jslaset = 2 
    1228             CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) ) 
    1229  
    1230          ENDIF 
    1231  
    1232          IF ( ln_slafb ) THEN 
    1233              
    1234             jfbini = jslaset + 1 
    1235  
    1236             DO jslaset = jfbini, nslasets 
    1237                 
    1238                jset = jslaset - jfbini + 1 
    1239                WRITE(cdtmp,'(A,I2.2)')'slafb_',jset 
    1240                CALL obs_wri_sla( cdtmp, sladata(jslaset) ) 
    1241  
    1242             END DO 
    1243  
    1244          ENDIF 
    1245  
    1246       ENDIF 
    1247  
    1248       !  - Sea surface temperature 
    1249       IF ( ln_sst ) THEN 
    1250  
    1251          ! Copy data from sstdatqc to sstdata structures 
    1252          DO jsstset = 1, nsstsets 
    1253       
    1254               CALL obs_surf_decompress( sstdatqc(jsstset), & 
    1255                  &                    sstdata(jsstset), .TRUE., numout ) 
    1256  
    1257          END DO 
    1258  
    1259          jsstset = 0  
    1260  
    1261          ! Write the AVISO SST data 
    1262  
    1263          IF ( ln_reysst ) THEN 
    1264              
    1265             jsstset = jsstset + 1 
    1266             CALL obs_wri_sst( 'reynolds', sstdata(jsstset) ) 
    1267  
    1268          ENDIF 
    1269  
    1270          IF ( ln_ghrsst ) THEN 
    1271              
    1272             jsstset = jsstset + 1 
    1273             CALL obs_wri_sst( 'ghr', sstdata(jsstset) ) 
    1274  
    1275          ENDIF 
    1276  
    1277          IF ( ln_sstfb ) THEN 
    1278              
    1279             jfbini = jsstset + 1 
    1280  
    1281             DO jsstset = jfbini, nsstsets 
    1282                 
    1283                jset = jsstset - jfbini + 1 
    1284                WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset 
    1285                CALL obs_wri_sst( cdtmp, sstdata(jsstset) ) 
    1286  
    1287             END DO 
    1288  
    1289          ENDIF 
    1290  
    1291       ENDIF 
    1292  
    1293       !  - Sea surface salinity 
    1294       IF ( ln_sss ) THEN 
    1295          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    1296       ENDIF 
    1297  
    1298       !  - Sea Ice Concentration 
    1299       IF ( ln_seaice ) THEN 
    1300  
    1301          ! Copy data from seaicedatqc to seaicedata structures 
    1302          DO jseaiceset = 1, nseaicesets 
    1303  
    1304               CALL obs_surf_decompress( seaicedatqc(jseaiceset), & 
    1305                  &                    seaicedata(jseaiceset), .TRUE., numout ) 
    1306  
    1307          END DO 
    1308  
    1309          ! Write the Sea Ice data 
    1310          DO jseaiceset = 1, nseaicesets 
    1311        
    1312             WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset 
    1313             CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) ) 
    1314  
    1315          END DO 
    1316  
    1317       ENDIF 
    1318        
    1319       ! Velocity data 
    1320       IF( ln_vel3d ) THEN 
    1321  
    1322          ! Copy data from veldatqc to velodata structures 
    1323          DO jveloset = 1, nvelosets 
    1324  
    1325             CALL obs_prof_decompress( veldatqc(jveloset), & 
    1326                  &                    velodata(jveloset), .TRUE., numout ) 
    1327  
    1328          END DO 
    1329  
    1330          ! Write the profiles. 
    1331  
    1332          jveloset = 0 
    1333  
    1334          ! Daily averaged data 
    1335  
    1336          IF ( ln_velavcur ) THEN 
    1337              
    1338             jveloset = jveloset + 1 
    1339  
    1340             CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint ) 
    1341  
    1342          ENDIF 
    1343  
    1344          ! High frequency data 
    1345  
    1346          IF ( ln_velhrcur ) THEN 
    1347              
    1348             jveloset = jveloset + 1 
    1349  
    1350             CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint ) 
    1351  
    1352          ENDIF 
    1353  
    1354          ! Daily averaged data 
    1355  
    1356          IF ( ln_velavadcp ) THEN 
    1357              
    1358             jveloset = jveloset + 1 
    1359  
    1360             CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint ) 
    1361  
    1362          ENDIF 
    1363  
    1364          ! High frequency data 
    1365  
    1366          IF ( ln_velhradcp ) THEN 
    1367              
    1368             jveloset = jveloset + 1 
    1369              
    1370             CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint ) 
    1371                 
    1372          ENDIF 
    1373  
    1374          ! Feedback velocity data 
    1375  
    1376          IF ( ln_velfb ) THEN 
    1377  
    1378             jfbini = jveloset + 1 
    1379  
    1380             DO jveloset = jfbini, nvelosets 
    1381                 
    1382                jset = jveloset - jfbini + 1 
    1383                WRITE(cdtmp,'(A,I2.2)')'velfb_',jset 
    1384                CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint ) 
    1385  
    1386             END DO 
    1387  
    1388          ENDIF 
    1389           
    1390       ENDIF 
     592          
     593      ENDIF 
     594 
    1391595 
    1392596   END SUBROUTINE dia_obs_wri 
     
    1409613 
    1410614      !! diaobs deallocation 
    1411       IF ( nprofsets > 0 ) THEN 
    1412           DEALLOCATE(ld_enact, & 
    1413                   &  profdata, & 
    1414                   &  prodatqc) 
    1415       END IF 
    1416       IF ( ln_sla ) THEN 
    1417           DEALLOCATE(sladata, & 
    1418                   &  sladatqc) 
    1419       END IF 
    1420       IF ( ln_seaice ) THEN 
    1421           DEALLOCATE(sladata, & 
    1422                   &  sladatqc) 
    1423       END IF 
    1424       IF ( ln_sst ) THEN 
    1425           DEALLOCATE(sstdata, & 
    1426                   &  sstdatqc) 
    1427       END IF 
    1428       IF ( ln_vel3d ) THEN 
    1429           DEALLOCATE(ld_velav, & 
    1430                   &  velodata, & 
    1431                   &  veldatqc) 
    1432       END IF 
     615      IF ( numproftypes > 0 ) DEALLOCATE(profdata, profdataqc, nvarsprof, nextrprof) 
     616      IF ( numsurftypes > 0 ) DEALLOCATE(surfdata, surfdataqc, nvarssurf, nextrsurf) 
     617       
    1433618   END SUBROUTINE dia_obs_dealloc 
    1434619 
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/julian.F90

    r2715 r5659  
    2424      &   greg2jul            ! Convert date to relative time  
    2525   
     26   !! $Id$ 
    2627CONTAINS 
    2728  
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r4245 r5659  
    77 
    88   !!---------------------------------------------------------------------- 
    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 
    2211   !!---------------------------------------------------------------------- 
    2312 
     
    4635   PRIVATE 
    4736 
    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 
     37   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile observations 
     38      &   obs_surf_opt     ! Compute the model counterpart of surface observations 
    5439 
    5540   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 
     
    6348CONTAINS 
    6449 
    65    SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    66       &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, & 
    67       &                    kdailyavtypes ) 
     50   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
     51      &                     pvar1, pvar2, pgdept, pmask1, k1dint, k2dint, & 
     52      &                     kdailyavtypes ) 
    6853      !!----------------------------------------------------------------------- 
    6954      !! 
     
    7863      !! 
    7964      !!    First, a vertical profile of horizontally interpolated model 
    80       !!    now temperatures is computed at the obs (lon, lat) point. 
     65      !!    now values is computed at the obs (lon, lat) point. 
    8166      !!    Several horizontal interpolation schemes are available: 
    8267      !!        - distance-weighted (great circle) (k2dint = 0) 
     
    8671      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    8772      !! 
    88       !!    Next, the vertical temperature profile is interpolated to the 
     73      !!    Next, the vertical profile is interpolated to the 
    8974      !!    data depth points. Two vertical interpolation schemes are 
    9075      !!    available: 
     
    9681      !!    routine. 
    9782      !! 
    98       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is 
     83      !!    If the logical is switched on, the model equivalent is 
    9984      !!    a daily mean model temperature field. So, we first compute 
    10085      !!    the mean, then interpolate only at the end of the day. 
    10186      !! 
    102       !!    Note: the in situ temperature observations must be converted 
     87      !!    Note: in situ temperature observations must be converted 
    10388      !!    to potential temperature (the model variable) prior to 
    10489      !!    assimilation.  
    105       !!?????????????????????????????????????????????????????????????? 
    106       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 
    107       !!?????????????????????????????????????????????????????????????? 
    10890      !! 
    10991      !! ** Action  : 
     
    11597      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity 
    11698      !!      ! 07-03 (K. Mogensen) General handling of profiles 
     99      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
    117100      !!----------------------------------------------------------------------- 
    118101   
     
    134117      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
    135118      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    136          & ptn,    &    ! Model temperature field 
    137          & psn,    &    ! Model salinity field 
    138          & ptmask       ! Land-sea mask 
     119         & pvar1,    &    ! Model field 1 
     120         & pvar2,    &    ! Model field 2 
     121         & pmask1,   &    ! Land-sea mask 1 
     122         & pmask2         ! Land-sea mask 2 
    139123      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    140124         & pgdept       ! Model array of depth levels 
     
    158142      REAL(KIND=wp) :: zdaystp 
    159143      REAL(KIND=wp), DIMENSION(kpk) :: & 
     144         & zobsmask1, & 
     145         & zobsmask2, & 
    160146         & zobsmask, & 
    161147         & zobsk,    & 
    162148         & zobs2k 
    163149      REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
    164          & zweig 
     150         & zweig1, & 
     151         & zweig2 
    165152      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    166          & zmask, & 
    167          & zintt, & 
    168          & zints, & 
    169          & zinmt, & 
    170          & zinms 
     153         & zmask1, & 
     154         & zmask2, & 
     155         & zint1, & 
     156         & zint2, & 
     157         & zinm1, & 
     158         & zinm2 
    171159      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    172          & zglam, & 
    173          & zgphi 
     160         & zglam1, & 
     161         & zglam2, & 
     162         & zgphi1, & 
     163         & zgphi2 
    174164      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    175          & igrdi, & 
    176          & igrdj 
     165         & igrdi1, & 
     166         & igrdi2, & 
     167         & igrdj1, & 
     168         & igrdj2 
    177169 
    178170      !------------------------------------------------------------------------ 
     
    209201         DO jj = 1, jpj 
    210202            DO ji = 1, jpi 
    211                ! Increment the temperature field for computing daily mean 
     203               ! Increment field 1 for computing daily mean 
    212204               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 
     205                  &                        + pvar1(ji,jj,jk) 
     206               ! Increment field 2 for computing daily mean 
    215207               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    216                   &                        + psn(ji,jj,jk) 
     208                  &                        + pvar2(ji,jj,jk) 
    217209            END DO 
    218210         END DO 
     
    236228      ! Get the data for interpolation 
    237229      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)   & 
     230         & igrdi1(2,2,ipro),      & 
     231         & igrdi2(2,2,ipro),      & 
     232         & igrdj1(2,2,ipro),      & 
     233         & igrdj2(2,2,ipro),      & 
     234         & zglam1(2,2,ipro),      & 
     235         & zglam2(2,2,ipro),      & 
     236         & zgphi1(2,2,ipro),      & 
     237         & zgphi2(2,2,ipro),      & 
     238         & zmask1(2,2,kpk,ipro),  & 
     239         & zmask2(2,2,kpk,ipro),  & 
     240         & zint1(2,2,kpk,ipro),  & 
     241         & zint2(2,2,kpk,ipro)   & 
    245242         & ) 
    246243 
    247244      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    248245         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) 
     246         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 
     247         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 
     248         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 
     249         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 
     250         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 
     251         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 
     252         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
     253         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
     254         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
     255         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
     256         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
     257         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
     258         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
     259         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
     260         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
     261         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
    257262      END DO 
    258263 
    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 ) 
     264      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, glamt1, zglam1 ) 
     265      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, gphit1, zgphi1 ) 
     266      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 
     267      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
     268       
     269      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, glamt2, zglam2 ) 
     270      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, gphit2, zgphi2 ) 
     271      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask1, zmask2 ) 
     272      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    264273 
    265274      ! At the end of the day also get interpolated means 
     
    267276 
    268277         ALLOCATE( & 
    269             & zinmt(2,2,kpk,ipro),  & 
    270             & zinms(2,2,kpk,ipro)   & 
     278            & zinm1(2,2,kpk,ipro),  & 
     279            & zinm2(2,2,kpk,ipro)   & 
    271280            & ) 
    272281 
    273          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    274             &                  prodatqc%vdmean(:,:,:,1), zinmt ) 
    275          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    276             &                  prodatqc%vdmean(:,:,:,2), zinms ) 
     282         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, & 
     283            &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
     284         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, & 
     285            &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
    277286 
    278287      ENDIF 
     
    304313         ! Horizontal weights and vertical mask 
    305314 
    306          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
    307             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 
     315         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    308316 
    309317            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    310                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    311                &                   zmask(:,:,:,iobs), zweig, zobsmask ) 
    312  
     318               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
     319               &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     320 
     321         ENDIF 
     322 
     323         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     324 
     325            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     326               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
     327               &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     328  
    313329         ENDIF 
    314330 
     
    321337               IF ( idayend == 0 )  THEN 
    322338                   
    323                   ! Daily averaged moored buoy (MRB) data 
     339                  ! Daily averaged data 
    324340                   
    325341                  CALL obs_int_h2d( kpk, kpk,      & 
    326                      &              zweig, zinmt(:,:,:,iobs), zobsk ) 
     342                     &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    327343                   
    328344                   
     
    330346                
    331347                  CALL ctl_stop( ' A nonzero' //     & 
    332                      &           ' number of profile T BUOY data should' // & 
     348                     &           ' number of average profile data should' // & 
    333349                     &           ' only occur at the end of a given day' ) 
    334350 
     
    340356 
    341357               CALL obs_int_h2d( kpk, kpk,      & 
    342                   &              zweig, zintt(:,:,:,iobs), zobsk ) 
     358                  &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    343359 
    344360            ENDIF 
     
    377393               IF ( idayend == 0 )  THEN 
    378394 
    379                   ! Daily averaged moored buoy (MRB) data 
     395                  ! Daily averaged data 
    380396                   
    381397                  CALL obs_int_h2d( kpk, kpk,      & 
    382                      &              zweig, zinms(:,:,:,iobs), zobsk ) 
     398                     &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    383399                   
    384400               ELSE 
    385401 
    386402                  CALL ctl_stop( ' A nonzero' //     & 
    387                      &           ' number of profile S BUOY data should' // & 
     403                     &           ' number of average profile data should' // & 
    388404                     &           ' only occur at the end of a given day' ) 
    389405 
     
    395411 
    396412               CALL obs_int_h2d( kpk, kpk,      & 
    397                   &              zweig, zints(:,:,:,iobs), zobsk ) 
     413                  &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    398414 
    399415            ENDIF 
     
    429445      ! Deallocate the data for interpolation 
    430446      DEALLOCATE( & 
    431          & igrdi, & 
    432          & igrdj, & 
    433          & zglam, & 
    434          & zgphi, & 
    435          & zmask, & 
    436          & zintt, & 
    437          & zints  & 
     447         & igrdi1, & 
     448         & igrdi2, & 
     449         & igrdj1, & 
     450         & igrdj2, & 
     451         & zglam1, & 
     452         & zglam2, & 
     453         & zgphi1, & 
     454         & zgphi2, & 
     455         & zmask1, & 
     456         & zmask2, & 
     457         & zint1,  & 
     458         & zint2   & 
    438459         & ) 
    439460      ! At the end of the day also get interpolated means 
    440461      IF ( idayend == 0 ) THEN 
    441462         DEALLOCATE( & 
    442             & zinmt,  & 
    443             & zinms   & 
     463            & zinm1,  & 
     464            & zinm2   & 
    444465            & ) 
    445466      ENDIF 
     
    447468      prodatqc%nprofup = prodatqc%nprofup + ipro  
    448469       
    449    END SUBROUTINE obs_pro_opt 
    450  
    451    SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 
    452       &                    psshn, psshmask, k2dint ) 
     470   END SUBROUTINE obs_prof_opt 
     471 
     472   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, kit000, & 
     473      &                    psurf, psurfmask, k2dint, ld_nightav ) 
    453474      !!----------------------------------------------------------------------- 
    454475      !! 
    455       !!                     ***  ROUTINE obs_sla_opt  *** 
    456       !! 
    457       !! ** Purpose : Compute the model counterpart of sea level anomaly 
     476      !!                     ***  ROUTINE obs_surf_opt  *** 
     477      !! 
     478      !! ** Purpose : Compute the model counterpart of surface 
    458479      !!              data by interpolating from the model grid to the  
    459480      !!              observation point. 
     
    462483      !!              the model values at the corners of the surrounding grid box. 
    463484      !! 
    464       !!    The now model SSH is first computed at the obs (lon, lat) point. 
     485      !!    The new model value is first computed at the obs (lon, lat) point. 
    465486      !! 
    466487      !!    Several horizontal interpolation schemes are available: 
     
    471492      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    472493      !!   
    473       !!    The sea level anomaly at the observation points is then computed  
    474       !!    by removing a mean dynamic topography (defined at the obs. point). 
    475494      !! 
    476495      !! ** Action  : 
     
    478497      !! History : 
    479498      !!      ! 07-03 (A. Weaver) 
     499      !!      ! 15-02 (M. Martin) Combined routine for surface types 
    480500      !!----------------------------------------------------------------------- 
    481501   
     
    486506 
    487507      !! * Arguments 
    488       TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening 
     508      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc     ! Subset of surface data not failing screening 
    489509      INTEGER, INTENT(IN) :: kt      ! Time step 
    490510      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
     
    494514      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    495515      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    496          & psshn,  &    ! Model SSH field 
    497          & psshmask     ! Land-sea mask 
     516         & psurf,  &    ! Model surface field 
     517         & psurfmask    ! Land-sea mask 
    498518          
    499519      !! * Local declarations 
     
    502522      INTEGER :: jobs 
    503523      INTEGER :: inrc 
    504       INTEGER :: isla 
     524      INTEGER :: isurf 
    505525      INTEGER :: iobs 
    506526      REAL(KIND=wp) :: zlam 
     
    511531      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    512532         & zmask, & 
    513          & zsshl, & 
     533         & zsurf, & 
    514534         & zglam, & 
    515535         & zgphi 
     
    523543      ! ... Record and data counters 
    524544      inrc = kt - kit000 + 2 
    525       isla = sladatqc%nsstp(inrc) 
     545      isurf = surfdataqc%nsstp(inrc) 
     546       
     547      IF ( ld_nightav ) THEN 
     548 
     549      ! Initialize array for night mean 
     550         IF ( kt .EQ. 0 ) THEN 
     551            ALLOCATE ( icount_night(kpi,kpj) ) 
     552            ALLOCATE ( imask_night(kpi,kpj) ) 
     553            ALLOCATE ( zintmp(kpi,kpj) ) 
     554            ALLOCATE ( zouttmp(kpi,kpj) ) 
     555            ALLOCATE ( zmeanday(kpi,kpj) ) 
     556            nday_qsr = -1   ! initialisation flag for nbc_dcy 
     557         ENDIF 
     558 
     559         ! Initialize daily mean for first timestep 
     560         idayend = MOD( kt - kit000 + 1, kdaystp ) 
     561 
     562         ! Added kt == 0 test to catch restart case  
     563         IF ( idayend == 1 .OR. kt == 0) THEN 
     564            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 
     565            DO jj = 1, jpj 
     566               DO ji = 1, jpi 
     567                  surfdataqc%vdmean(ji,jj) = 0.0 
     568                  zmeanday(ji,jj) = 0.0 
     569                  icount_night(ji,jj) = 0 
     570               END DO 
     571            END DO 
     572         ENDIF 
     573 
     574         zintmp(:,:) = 0.0 
     575         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
     576         imask_night(:,:) = INT( zouttmp(:,:) ) 
     577 
     578         DO jj = 1, jpj 
     579            DO ji = 1, jpi 
     580               ! Increment the temperature field for computing night mean and counter 
     581               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  & 
     582                      &                        + psurf(ji,jj)*imask_night(ji,jj) 
     583               zmeanday(ji,jj)        = zmeanday(ji,jj) + psurf(ji,jj) 
     584               icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 
     585            END DO 
     586         END DO 
     587    
     588         ! Compute the daily mean at the end of day 
     589 
     590         zdaystp = 1.0 / REAL( kdaystp ) 
     591 
     592         IF ( idayend == 0 ) THEN  
     593            DO jj = 1, jpj 
     594               DO ji = 1, jpi 
     595                  ! Test if "no night" point 
     596                  IF ( icount_night(ji,jj) .NE. 0 ) THEN 
     597                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     598                       &                        / icount_night(ji,jj)  
     599                  ELSE 
     600                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
     601                  ENDIF 
     602               END DO 
     603            END DO 
     604         ENDIF 
     605 
     606      ENDIF 
    526607 
    527608      ! Get the data for interpolation 
    528609 
    529610      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)  & 
     611         & igrdi(2,2,isurf), & 
     612         & igrdj(2,2,isurf), & 
     613         & zglam(2,2,isurf), & 
     614         & zgphi(2,2,isurf), & 
     615         & zmask(2,2,isurf), & 
     616         & zsurf(2,2,isurf)  & 
    536617         & ) 
    537618       
    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) 
     619      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     620         iobs = jobs - surfdataqc%nsurfup 
     621         igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
     622         igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
     623         igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
     624         igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
     625         igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
     626         igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
     627         igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
     628         igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
    548629      END DO 
    549630 
    550       CALL obs_int_comm_2d( 2, 2, isla, & 
     631      CALL obs_int_comm_2d( 2, 2, isurf, & 
    551632         &                  igrdi, igrdj, glamt, zglam ) 
    552       CALL obs_int_comm_2d( 2, 2, isla, & 
     633      CALL obs_int_comm_2d( 2, 2, isurf, & 
    553634         &                  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 ) 
     635      CALL obs_int_comm_2d( 2, 2, isurf, & 
     636         &                  igrdi, igrdj, psurfmask, zmask ) 
     637      CALL obs_int_comm_2d( 2, 2, isurf, & 
     638         &                  igrdi, igrdj, psurf, zsurf ) 
     639 
     640      ! At the end of the day get interpolated means 
     641      IF ( idayend == 0 .AND. ld_nightav ) THEN 
     642 
     643         ALLOCATE( & 
     644            & zsurfm(2,2,isurf)  & 
     645            & ) 
     646 
     647         CALL obs_int_comm_2d( 2, 2, isurf, igrdi, igrdj, & 
     648            &               surfdataqc%vdmean(:,:), zsurfm ) 
     649 
     650      ENDIF 
    558651 
    559652      ! 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 
     653      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     654 
     655         iobs = jobs - surfdataqc%nsurfup 
     656 
     657         IF ( kt /= surfdataqc%mstp(jobs) ) THEN 
    566658             
    567659            IF(lwp) THEN 
     
    574666               WRITE(numout,*) ' Record  = ', jobs,                & 
    575667                  &            ' kt      = ', kt,                  & 
    576                   &            ' mstp    = ', sladatqc%mstp(jobs), & 
    577                   &            ' ntyp    = ', sladatqc%ntyp(jobs) 
     668                  &            ' mstp    = ', surfdataqc%mstp(jobs), & 
     669                  &            ' ntyp    = ', surfdataqc%ntyp(jobs) 
    578670            ENDIF 
    579671            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 
     
    581673         ENDIF 
    582674          
    583          zlam = sladatqc%rlam(jobs) 
    584          zphi = sladatqc%rphi(jobs) 
    585  
    586          ! Get weights to interpolate the model SSH to the observation point 
     675         zlam = surfdataqc%rlam(jobs) 
     676         zphi = surfdataqc%rphi(jobs) 
     677 
     678         ! Get weights to interpolate the model value to the observation point 
    587679         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    588680            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     
    590682          
    591683 
    592          ! Interpolate the model SSH to the observation point 
    593          CALL obs_int_h2d( 1, 1,      & 
    594             &              zweig, zsshl(:,:,iobs),  zext ) 
     684         ! Interpolate the model field to the observation point 
     685         IF ( ld_nightav ) THEN 
     686 
     687           IF ( idayend == 0 )  THEN 
     688               ! Daily averaged data 
     689               CALL obs_int_h2d( 1, 1,      &  
     690                     &              zweig, zsurfm(:,:,iobs), zext ) 
     691            ELSE  
     692               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     & 
     693                     &           ' number of night data should' // & 
     694                     &           ' only occur at the end of a given day' ) 
     695            ENDIF 
     696 
     697         ELSE 
     698 
     699            CALL obs_int_h2d( 1, 1,      & 
     700               &              zweig, zsurf(:,:,iobs),  zext ) 
     701 
     702         ENDIF 
    595703          
    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) 
     704         IF ( surfdataqc%nextr == 2 ) THEN 
     705            ! ... Remove the MDT from the SSH at the observation point to get the SLA 
     706            surfdataqc%rext(jobs,1) = zext(1) 
     707            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
     708         ELSE 
     709            surfdataqc%rmod(jobs,1) = zext(1) 
     710         ENDIF 
    599711 
    600712      END DO 
     
    607719         & zgphi, & 
    608720         & zmask, & 
    609          & zsshl  & 
    610          & ) 
    611  
    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  & 
     721         & zsurf  & 
    875722         & ) 
    876723 
     
    878725      IF ( idayend == 0 .AND. ld_nightav ) THEN 
    879726         DEALLOCATE( & 
    880             & zsstm  & 
     727            & zsurfm  & 
    881728            & ) 
    882729      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 
     730 
     731      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     732 
     733   END SUBROUTINE obs_surf_opt 
    1440734 
    1441735END MODULE obs_oper 
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r4292 r5659  
    77 
    88   !!--------------------------------------------------------------------- 
    9    !!   obs_pre_pro  : First level check and screening of T/S profiles 
    10    !!   obs_pre_sla  : First level check and screening of SLA observations 
    11    !!   obs_pre_sst  : First level check and screening of SLA observations 
    12    !!   obs_pre_seaice : First level check and screening of sea ice observations 
    13    !!   obs_pre_vel  : First level check and screening of velocity obs. 
    14    !!   obs_scr      : Basic screening of the observations 
    15    !!   obs_coo_tim  : Compute number of time steps to the observation time 
    16    !!   obs_sor      : Sort the observation arrays 
     9   !!   obs_pre_prof  : First level check and screening of profile observations 
     10   !!   obs_pre_surf  : First level check and screening of surface observations 
     11   !!   obs_scr       : Basic screening of the observations 
     12   !!   obs_coo_tim   : Compute number of time steps to the observation time 
     13   !!   obs_sor       : Sort the observation arrays 
    1714   !!--------------------------------------------------------------------- 
    1815   !! * Modules used 
     
    3633 
    3734   PUBLIC & 
    38       & obs_pre_pro, &    ! First level check and screening of profiles 
    39       & obs_pre_sla, &    ! First level check and screening of SLA data 
    40       & obs_pre_sst, &    ! First level check and screening of SLA data 
    41       & obs_pre_seaice, & ! First level check and screening of sea ice data 
    42       & obs_pre_vel, &     ! First level check and screening of velocity profiles 
    43       & calc_month_len     ! Calculate the number of days in the months of a year   
     35      & obs_pre_prof, &    ! First level check and screening of profile obs 
     36      & obs_pre_surf, &    ! First level check and screening of surface obs 
     37      & calc_month_len     ! Calculate the number of days in the months of a year 
    4438 
    4539   !!---------------------------------------------------------------------- 
     
    5145CONTAINS 
    5246 
    53    SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, & 
    54       &                    kdailyavtypes ) 
    55       !!---------------------------------------------------------------------- 
    56       !!                    ***  ROUTINE obs_pre_pro  *** 
    57       !! 
    58       !! ** Purpose : First level check and screening of T and S profiles 
    59       !! 
    60       !! ** Method  : First level check and screening of T and S profiles 
    61       !! 
    62       !! ** Action  :  
    63       !! 
    64       !! References : 
    65       !!    
    66       !! History : 
    67       !!        !  2007-01  (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d  
    68       !!        !  2007-03  (K. Mogensen) General handling of profiles 
    69       !!        !  2007-06  (K. Mogensen et al) Reject obs. near land. 
    70       !!---------------------------------------------------------------------- 
    71       !! * Modules used 
    72       USE domstp              ! Domain: set the time-step 
    73       USE par_oce             ! Ocean parameters 
    74       USE dom_oce, ONLY : &   ! Geographical information 
    75          & glamt,   & 
    76          & gphit,   & 
    77          & gdept_1d,& 
    78          & tmask,   & 
    79          & nproc 
    80       !! * Arguments 
    81       TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data 
    82       TYPE(obs_prof), INTENT(INOUT) :: prodatqc     ! Subset of profile data not failing screening 
    83       LOGICAL, INTENT(IN) :: ld_t3d         ! Switch for temperature 
    84       LOGICAL, INTENT(IN) :: ld_s3d         ! Switch for salinity 
    85       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    86       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    87          & kdailyavtypes! Types for daily averages 
    88       !! * Local declarations    
    89       INTEGER :: iyea0         ! Initial date 
    90       INTEGER :: imon0         !  - (year, month, day, hour, minute) 
    91       INTEGER :: iday0    
    92       INTEGER :: ihou0 
    93       INTEGER :: imin0 
    94       INTEGER :: icycle        ! Current assimilation cycle 
    95                                ! Counters for observations that 
    96       INTEGER :: iotdobs       !  - outside time domain 
    97       INTEGER :: iosdtobs      !  - outside space domain (temperature) 
    98       INTEGER :: iosdsobs      !  - outside space domain (salinity) 
    99       INTEGER :: ilantobs      !  - within a model land cell (temperature) 
    100       INTEGER :: ilansobs      !  - within a model land cell (salinity) 
    101       INTEGER :: inlatobs      !  - close to land (temperature) 
    102       INTEGER :: inlasobs      !  - close to land (salinity) 
    103       INTEGER :: igrdobs       !  - fail the grid search 
    104                                ! Global counters for observations that 
    105       INTEGER :: iotdobsmpp    !  - outside time domain 
    106       INTEGER :: iosdtobsmpp   !  - outside space domain (temperature) 
    107       INTEGER :: iosdsobsmpp   !  - outside space domain (salinity) 
    108       INTEGER :: ilantobsmpp   !  - within a model land cell (temperature) 
    109       INTEGER :: ilansobsmpp   !  - within a model land cell (salinity) 
    110       INTEGER :: inlatobsmpp   !  - close to land (temperature) 
    111       INTEGER :: inlasobsmpp   !  - close to land (salinity) 
    112       INTEGER :: igrdobsmpp    !  - fail the grid search 
    113       TYPE(obs_prof_valid) ::  llvalid     ! Profile selection  
    114       TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 
    115          & llvvalid            ! T,S selection  
    116       INTEGER :: jvar          ! Variable loop variable 
    117       INTEGER :: jobs          ! Obs. loop variable 
    118       INTEGER :: jstp          ! Time loop variable 
    119       INTEGER :: inrc          ! Time index variable 
    120        
    121       IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...' 
    122  
    123       ! Initial date initialization (year, month, day, hour, minute) 
    124       iyea0 =   ndate0 / 10000 
    125       imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
    126       iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
    127       ihou0 = 0 
    128       imin0 = 0 
    129  
    130       icycle = no     ! Assimilation cycle 
    131  
    132       ! Diagnotics counters for various failures. 
    133  
    134       iotdobs  = 0 
    135       igrdobs  = 0 
    136       iosdtobs = 0 
    137       iosdsobs = 0 
    138       ilantobs = 0 
    139       ilansobs = 0 
    140       inlatobs = 0 
    141       inlasobs = 0 
    142  
    143       ! ----------------------------------------------------------------------- 
    144       ! Find time coordinate for profiles 
    145       ! ----------------------------------------------------------------------- 
    146  
    147       IF ( PRESENT(kdailyavtypes) ) THEN 
    148          CALL obs_coo_tim_prof( icycle, & 
    149             &                iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    150             &                profdata%nprof,   profdata%nyea, profdata%nmon, & 
    151             &                profdata%nday,    profdata%nhou, profdata%nmin, & 
    152             &                profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    153             &                iotdobs, kdailyavtypes = kdailyavtypes        ) 
    154       ELSE 
    155          CALL obs_coo_tim_prof( icycle, & 
    156             &                iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    157             &                profdata%nprof,   profdata%nyea, profdata%nmon, & 
    158             &                profdata%nday,    profdata%nhou, profdata%nmin, & 
    159             &                profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    160             &                iotdobs ) 
    161       ENDIF 
    162       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
    163        
    164       ! ----------------------------------------------------------------------- 
    165       ! Check for profiles failing the grid search 
    166       ! ----------------------------------------------------------------------- 
    167  
    168       CALL obs_coo_grd( profdata%nprof,   profdata%mi, profdata%mj, & 
    169          &              profdata%nqc,     igrdobs                         ) 
    170  
    171       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    172  
    173       ! ----------------------------------------------------------------------- 
    174       ! Reject all observations for profiles with nqc > 10 
    175       ! ----------------------------------------------------------------------- 
    176  
    177       CALL obs_pro_rej( profdata ) 
    178  
    179       ! ----------------------------------------------------------------------- 
    180       ! Check for land points. This includes points below the model 
    181       ! bathymetry so this is done for every point in the profile 
    182       ! ----------------------------------------------------------------------- 
    183  
    184       ! Temperature 
    185  
    186       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   & 
    187          &                 profdata%npvsta(:,1),  profdata%npvend(:,1), & 
    188          &                 jpi,                   jpj,                  & 
    189          &                 jpk,                                         & 
    190          &                 profdata%mi,           profdata%mj,          &  
    191          &                 profdata%var(1)%mvk,                         & 
    192          &                 profdata%rlam,         profdata%rphi,        & 
    193          &                 profdata%var(1)%vdep,                        & 
    194          &                 glamt,                 gphit,                & 
    195          &                 gdept_1d,              tmask,                & 
    196          &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    197          &                 iosdtobs,              ilantobs,             & 
    198          &                 inlatobs,              ld_nea                ) 
    199  
    200       CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp ) 
    201       CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp ) 
    202       CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp ) 
    203  
    204       ! Salinity 
    205  
    206       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   & 
    207          &                 profdata%npvsta(:,2),  profdata%npvend(:,2), & 
    208          &                 jpi,                   jpj,                  & 
    209          &                 jpk,                                         & 
    210          &                 profdata%mi,           profdata%mj,          &  
    211          &                 profdata%var(2)%mvk,                         & 
    212          &                 profdata%rlam,         profdata%rphi,        & 
    213          &                 profdata%var(2)%vdep,                        & 
    214          &                 glamt,                 gphit,                & 
    215          &                 gdept_1d,              tmask,                & 
    216          &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    217          &                 iosdsobs,              ilansobs,             & 
    218          &                 inlasobs,              ld_nea                ) 
    219  
    220       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    221       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    222       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
    223  
    224       ! ----------------------------------------------------------------------- 
    225       ! Copy useful data from the profdata data structure to 
    226       ! the prodatqc data structure  
    227       ! ----------------------------------------------------------------------- 
    228  
    229       ! Allocate the selection arrays 
    230  
    231       ALLOCATE( llvalid%luse(profdata%nprof) ) 
    232       DO jvar = 1,profdata%nvar 
    233          ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) ) 
    234       END DO 
    235  
    236       ! We want all data which has qc flags <= 10 
    237  
    238       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
    239       DO jvar = 1,profdata%nvar 
    240          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
    241       END DO 
    242  
    243       ! The actual copying 
    244  
    245       CALL obs_prof_compress( profdata,     prodatqc,       .TRUE.,  numout, & 
    246          &                    lvalid=llvalid, lvvalid=llvvalid ) 
    247  
    248       ! Dellocate the selection arrays 
    249       DEALLOCATE( llvalid%luse ) 
    250       DO jvar = 1,profdata%nvar 
    251          DEALLOCATE( llvvalid(jvar)%luse ) 
    252       END DO 
    253  
    254       ! ----------------------------------------------------------------------- 
    255       ! Print information about what observations are left after qc 
    256       ! ----------------------------------------------------------------------- 
    257  
    258       ! Update the total observation counter array 
    259        
    260       IF(lwp) THEN 
    261          WRITE(numout,*) 
    262          WRITE(numout,*) 'obs_pre_pro :' 
    263          WRITE(numout,*) '~~~~~~~~~~~' 
    264          WRITE(numout,*) 
    265          WRITE(numout,*) ' Profiles outside time domain                = ', & 
    266             &            iotdobsmpp 
    267          WRITE(numout,*) ' Remaining profiles that failed grid search  = ', & 
    268             &            igrdobsmpp 
    269          WRITE(numout,*) ' Remaining T data outside space domain       = ', & 
    270             &            iosdtobsmpp 
    271          WRITE(numout,*) ' Remaining T data at land points             = ', & 
    272             &            ilantobsmpp 
    273          IF (ld_nea) THEN 
    274             WRITE(numout,*) ' Remaining T data near land points (removed) = ',& 
    275                &            inlatobsmpp 
    276          ELSE 
    277             WRITE(numout,*) ' Remaining T data near land points (kept)    = ',& 
    278                &            inlatobsmpp 
    279          ENDIF 
    280          WRITE(numout,*) ' T data accepted                             = ', & 
    281             &            prodatqc%nvprotmpp(1) 
    282          WRITE(numout,*) ' Remaining S data outside space domain       = ', & 
    283             &            iosdsobsmpp 
    284          WRITE(numout,*) ' Remaining S data at land points             = ', & 
    285             &            ilansobsmpp 
    286          IF (ld_nea) THEN 
    287             WRITE(numout,*) ' Remaining S data near land points (removed) = ',& 
    288                &            inlasobsmpp 
    289          ELSE 
    290             WRITE(numout,*) ' Remaining S data near land points (kept)    = ',& 
    291                &            inlasobsmpp 
    292          ENDIF 
    293          WRITE(numout,*) ' S data accepted                             = ', & 
    294             &            prodatqc%nvprotmpp(2) 
    295  
    296          WRITE(numout,*) 
    297          WRITE(numout,*) ' Number of observations per time step :' 
    298          WRITE(numout,*) 
    299          WRITE(numout,997) 
    300          WRITE(numout,998) 
    301       ENDIF 
    302        
    303       DO jobs = 1, prodatqc%nprof 
    304          inrc = prodatqc%mstp(jobs) + 2 - nit000 
    305          prodatqc%npstp(inrc)  = prodatqc%npstp(inrc) + 1 
    306          DO jvar = 1, prodatqc%nvar 
    307             IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN 
    308                prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + & 
    309                   &                      ( prodatqc%npvend(jobs,jvar) - & 
    310                   &                        prodatqc%npvsta(jobs,jvar) + 1 ) 
    311             ENDIF 
    312          END DO 
    313       END DO 
    314        
    315        
    316       CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, & 
    317          &                       nitend - nit000 + 2 ) 
    318       DO jvar = 1, prodatqc%nvar 
    319          CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), & 
    320             &                       prodatqc%nvstpmpp(:,jvar), & 
    321             &                       nitend - nit000 + 2 ) 
    322       END DO 
    323  
    324       IF ( lwp ) THEN 
    325          DO jstp = nit000 - 1, nitend 
    326             inrc = jstp - nit000 + 2 
    327             WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
    328                &                    prodatqc%nvstpmpp(inrc,1), & 
    329                &                    prodatqc%nvstpmpp(inrc,2) 
    330          END DO 
    331       ENDIF 
    332  
    333 997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity') 
    334 998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------') 
    335 999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 
    336        
    337    END SUBROUTINE obs_pre_pro 
    338  
    339    SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea ) 
     47   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
    34048      !!---------------------------------------------------------------------- 
    34149      !!                    ***  ROUTINE obs_pre_sla  *** 
    34250      !! 
    343       !! ** Purpose : First level check and screening of SLA observations 
    344       !! 
    345       !! ** Method  : First level check and screening of SLA observations 
     51      !! ** Purpose : First level check and screening of surface observations 
     52      !! 
     53      !! ** Method  : First level check and screening of surface observations 
    34654      !! 
    34755      !! ** Action  :  
     
    35260      !!        !  2007-03  (A. Weaver, K. Mogensen) Original 
    35361      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land. 
     62      !!        !  2015-02  (M. Martin) Combined routine for surface types. 
    35463      !!---------------------------------------------------------------------- 
    35564      !! * Modules used 
     
    36271         & nproc 
    36372      !! * Arguments 
    364       TYPE(obs_surf), INTENT(INOUT) :: sladata    ! Full set of SLA data 
    365       TYPE(obs_surf), INTENT(INOUT) :: sladatqc   ! Subset of SLA data not failing screening 
    366       LOGICAL, INTENT(IN) :: ld_sla         ! Switch for SLA data 
     73      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of SLA data 
     74      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of SLA data not failing screening 
    36775      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    36876      !! * Local declarations 
     
    411119 
    412120      ! ----------------------------------------------------------------------- 
    413       ! Find time coordinate for SLA data 
     121      ! Find time coordinate for surface data 
    414122      ! ----------------------------------------------------------------------- 
    415123 
    416124      CALL obs_coo_tim( icycle, & 
    417125         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    418          &              sladata%nsurf,   sladata%nyea, sladata%nmon, & 
    419          &              sladata%nday,    sladata%nhou, sladata%nmin, & 
    420          &              sladata%nqc,     sladata%mstp, iotdobs        ) 
     126         &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
     127         &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
     128         &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
    421129 
    422130      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
    423131       
    424132      ! ----------------------------------------------------------------------- 
    425       ! Check for SLA data failing the grid search 
    426       ! ----------------------------------------------------------------------- 
    427  
    428       CALL obs_coo_grd( sladata%nsurf,   sladata%mi, sladata%mj, & 
    429          &              sladata%nqc,     igrdobs                         ) 
     133      ! Check for surface data failing the grid search 
     134      ! ----------------------------------------------------------------------- 
     135 
     136      CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, & 
     137         &              surfdata%nqc,     igrdobs                         ) 
    430138 
    431139      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     
    435143      ! ----------------------------------------------------------------------- 
    436144 
    437       CALL obs_coo_spc_2d( sladata%nsurf,              & 
     145      CALL obs_coo_spc_2d( surfdata%nsurf,              & 
    438146         &                 jpi,          jpj,          & 
    439          &                 sladata%mi,   sladata%mj,   &  
    440          &                 sladata%rlam, sladata%rphi, & 
     147         &                 surfdata%mi,   surfdata%mj,   &  
     148         &                 surfdata%rlam, surfdata%rphi, & 
    441149         &                 glamt,        gphit,        & 
    442          &                 tmask(:,:,1), sladata%nqc,  & 
     150         &                 tmask(:,:,1), surfdata%nqc,  & 
    443151         &                 iosdsobs,     ilansobs,     & 
    444152         &                 inlasobs,     ld_nea        ) 
     
    449157 
    450158      ! ----------------------------------------------------------------------- 
    451       ! Copy useful data from the sladata data structure to 
    452       ! the sladatqc data structure  
     159      ! Copy useful data from the surfdata data structure to 
     160      ! the surfdataqc data structure  
    453161      ! ----------------------------------------------------------------------- 
    454162 
    455163      ! Allocate the selection arrays 
    456164 
    457       ALLOCATE( llvalid(sladata%nsurf) ) 
     165      ALLOCATE( llvalid(surfdata%nsurf) ) 
    458166       
    459167      ! We want all data which has qc flags <= 10 
    460168 
    461       llvalid(:)  = ( sladata%nqc(:)  <= 10 ) 
     169      llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
    462170 
    463171      ! The actual copying 
    464172 
    465       CALL obs_surf_compress( sladata,     sladatqc,       .TRUE.,  numout, & 
     173      CALL obs_surf_compress( surfdata,     surfdataqc,       .TRUE.,  numout, & 
    466174         &                    lvalid=llvalid ) 
    467175 
     
    477185      IF(lwp) THEN 
    478186         WRITE(numout,*) 
    479          WRITE(numout,*) 'obs_pre_sla :' 
     187         WRITE(numout,*) 'obs_pre_surf :' 
    480188         WRITE(numout,*) '~~~~~~~~~~~' 
    481189         WRITE(numout,*) 
    482          WRITE(numout,*) ' SLA data outside time domain                  = ', & 
     190         WRITE(numout,*) ' Surface data outside time domain               = ', & 
    483191            &            iotdobsmpp 
    484          WRITE(numout,*) ' Remaining SLA data that failed grid search    = ', & 
     192         WRITE(numout,*) ' Remaining surface data that failed grid search    = ', & 
    485193            &            igrdobsmpp 
    486          WRITE(numout,*) ' Remaining SLA data outside space domain       = ', & 
     194         WRITE(numout,*) ' Remaining surface data outside space domain       = ', & 
    487195            &            iosdsobsmpp 
    488          WRITE(numout,*) ' Remaining SLA data at land points             = ', & 
     196         WRITE(numout,*) ' Remaining surface data at land points             = ', & 
    489197            &            ilansobsmpp 
    490198         IF (ld_nea) THEN 
    491             WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', & 
     199            WRITE(numout,*) ' Remaining surface data near land points (removed) = ', & 
    492200               &            inlasobsmpp 
    493201         ELSE 
    494             WRITE(numout,*) ' Remaining SLA data near land points (kept)    = ', & 
     202            WRITE(numout,*) ' Remaining surface data near land points (kept)    = ', & 
    495203               &            inlasobsmpp 
    496204         ENDIF 
    497          WRITE(numout,*) ' SLA data accepted                             = ', & 
    498             &            sladatqc%nsurfmpp 
     205         WRITE(numout,*) ' surface data accepted                             = ', & 
     206            &            surfdataqc%nsurfmpp 
    499207 
    500208         WRITE(numout,*) 
     
    505213      ENDIF 
    506214       
    507       DO jobs = 1, sladatqc%nsurf 
    508          inrc = sladatqc%mstp(jobs) + 2 - nit000 
    509          sladatqc%nsstp(inrc)  = sladatqc%nsstp(inrc) + 1 
     215      DO jobs = 1, surfdataqc%nsurf 
     216         inrc = surfdataqc%mstp(jobs) + 2 - nit000 
     217         surfdataqc%nsstp(inrc)  = surfdataqc%nsstp(inrc) + 1 
    510218      END DO 
    511219       
    512       CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, & 
     220      CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, & 
    513221         &                       nitend - nit000 + 2 ) 
    514222 
     
    516224         DO jstp = nit000 - 1, nitend 
    517225            inrc = jstp - nit000 + 2 
    518             WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc) 
     226            WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc) 
    519227         END DO 
    520228      ENDIF 
     
    5242321999  FORMAT(10X,I9,5X,I17) 
    525233 
    526    END SUBROUTINE obs_pre_sla 
    527  
    528    SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea ) 
    529       !!---------------------------------------------------------------------- 
    530       !!                    ***  ROUTINE obs_pre_sst  *** 
    531       !! 
    532       !! ** Purpose : First level check and screening of SST observations 
    533       !! 
    534       !! ** Method  : First level check and screening of SST observations 
    535       !! 
    536       !! ** Action  :  
    537       !! 
    538       !! References : 
    539       !!    
    540       !! History : 
    541       !!        !  2007-03  (S. Ricci) SST data preparation  
    542       !!---------------------------------------------------------------------- 
    543       !! * Modules used 
    544       USE domstp              ! Domain: set the time-step 
    545       USE par_oce             ! Ocean parameters 
    546       USE dom_oce, ONLY : &   ! Geographical information 
    547          & glamt,   & 
    548          & gphit,   & 
    549          & tmask,   & 
    550          & nproc 
    551       !! * Arguments 
    552       TYPE(obs_surf), INTENT(INOUT) :: sstdata     ! Full set of SST data 
    553       TYPE(obs_surf), INTENT(INOUT) :: sstdatqc    ! Subset of SST data not failing screening 
    554       LOGICAL :: ld_sst             ! Switch for SST data 
    555       LOGICAL :: ld_nea             ! Switch for rejecting observation near land 
    556       !! * Local declarations 
    557       INTEGER :: iyea0        ! Initial date 
    558       INTEGER :: imon0        !  - (year, month, day, hour, minute) 
    559       INTEGER :: iday0    
    560       INTEGER :: ihou0     
    561       INTEGER :: imin0 
    562       INTEGER :: icycle       ! Current assimilation cycle 
    563                               ! Counters for observations that 
    564       INTEGER :: iotdobs      !  - outside time domain 
    565       INTEGER :: iosdsobs     !  - outside space domain 
    566       INTEGER :: ilansobs     !  - within a model land cell 
    567       INTEGER :: inlasobs     !  - close to land 
    568       INTEGER :: igrdobs      !  - fail the grid search 
    569                               ! Global counters for observations that 
    570       INTEGER :: iotdobsmpp   !  - outside time domain 
    571       INTEGER :: iosdsobsmpp  !  - outside space domain 
    572       INTEGER :: ilansobsmpp  !  - within a model land cell 
    573       INTEGER :: inlasobsmpp  !  - close to land 
    574       INTEGER :: igrdobsmpp   !  - fail the grid search 
    575       LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    576          & llvalid            ! SST data selection 
    577       INTEGER :: jobs         ! Obs. loop variable 
    578       INTEGER :: jstp         ! Time loop variable 
    579       INTEGER :: inrc         ! Time index variable 
    580  
    581       IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...' 
    582  
    583       ! Initial date initialization (year, month, day, hour, minute) 
    584       iyea0 =   ndate0 / 10000 
    585       imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
    586       iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
    587       ihou0 = 0 
    588       imin0 = 0 
    589  
    590       icycle = no     ! Assimilation cycle 
    591  
    592       ! Diagnotics counters for various failures. 
    593  
    594       iotdobs  = 0 
    595       igrdobs  = 0 
    596       iosdsobs = 0 
    597       ilansobs = 0 
    598       inlasobs = 0 
    599  
    600       ! ----------------------------------------------------------------------- 
    601       ! Find time coordinate for SST data 
    602       ! ----------------------------------------------------------------------- 
    603  
    604       CALL obs_coo_tim( icycle, & 
    605          &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    606          &              sstdata%nsurf,   sstdata%nyea, sstdata%nmon, & 
    607          &              sstdata%nday,    sstdata%nhou, sstdata%nmin, & 
    608          &              sstdata%nqc,     sstdata%mstp, iotdobs        ) 
    609       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
    610       ! ----------------------------------------------------------------------- 
    611       ! Check for SST data failing the grid search 
    612       ! ----------------------------------------------------------------------- 
    613  
    614       CALL obs_coo_grd( sstdata%nsurf,   sstdata%mi, sstdata%mj, & 
    615          &              sstdata%nqc,     igrdobs                         ) 
    616       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    617  
    618       ! ----------------------------------------------------------------------- 
    619       ! Check for land points.  
    620       ! ----------------------------------------------------------------------- 
    621  
    622       CALL obs_coo_spc_2d( sstdata%nsurf,              & 
    623          &                 jpi,          jpj,          & 
    624          &                 sstdata%mi,   sstdata%mj,   &  
    625          &                 sstdata%rlam, sstdata%rphi, & 
    626          &                 glamt,        gphit,        & 
    627          &                 tmask(:,:,1), sstdata%nqc,  & 
    628          &                 iosdsobs,     ilansobs,     & 
    629          &                 inlasobs,     ld_nea        ) 
    630  
    631       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    632       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    633       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
    634  
    635       ! ----------------------------------------------------------------------- 
    636       ! Copy useful data from the sstdata data structure to 
    637       ! the sstdatqc data structure  
    638       ! ----------------------------------------------------------------------- 
    639  
    640       ! Allocate the selection arrays 
    641  
    642       ALLOCATE( llvalid(sstdata%nsurf) ) 
    643        
    644       ! We want all data which has qc flags <= 0 
    645  
    646       llvalid(:)  = ( sstdata%nqc(:)  <= 10 ) 
    647  
    648       ! The actual copying 
    649  
    650       CALL obs_surf_compress( sstdata,     sstdatqc,       .TRUE.,  numout, & 
    651          &                    lvalid=llvalid ) 
    652  
    653       ! Dellocate the selection arrays 
    654       DEALLOCATE( llvalid ) 
    655  
    656       ! ----------------------------------------------------------------------- 
    657       ! Print information about what observations are left after qc 
    658       ! ----------------------------------------------------------------------- 
    659  
    660       ! Update the total observation counter array 
    661        
    662       IF(lwp) THEN 
    663          WRITE(numout,*) 
    664          WRITE(numout,*) 'obs_pre_sst :' 
    665          WRITE(numout,*) '~~~~~~~~~~~' 
    666          WRITE(numout,*) 
    667          WRITE(numout,*) ' SST data outside time domain                  = ', & 
    668             &            iotdobsmpp 
    669          WRITE(numout,*) ' Remaining SST data that failed grid search    = ', & 
    670             &            igrdobsmpp 
    671          WRITE(numout,*) ' Remaining SST data outside space domain       = ', & 
    672             &            iosdsobsmpp 
    673          WRITE(numout,*) ' Remaining SST data at land points             = ', & 
    674             &            ilansobsmpp 
    675          IF (ld_nea) THEN 
    676             WRITE(numout,*) ' Remaining SST data near land points (removed) = ', & 
    677                &            inlasobsmpp 
    678          ELSE 
    679             WRITE(numout,*) ' Remaining SST data near land points (kept)    = ', & 
    680                &            inlasobsmpp 
    681          ENDIF 
    682          WRITE(numout,*) ' SST data accepted                             = ', & 
    683             &            sstdatqc%nsurfmpp 
    684  
    685          WRITE(numout,*) 
    686          WRITE(numout,*) ' Number of observations per time step :' 
    687          WRITE(numout,*) 
    688          WRITE(numout,1997) 
    689          WRITE(numout,1998) 
    690       ENDIF 
    691        
    692       DO jobs = 1, sstdatqc%nsurf 
    693          inrc = sstdatqc%mstp(jobs) + 2 - nit000 
    694          sstdatqc%nsstp(inrc)  = sstdatqc%nsstp(inrc) + 1 
    695       END DO 
    696        
    697       CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, & 
    698          &                       nitend - nit000 + 2 ) 
    699  
    700       IF ( lwp ) THEN 
    701          DO jstp = nit000 - 1, nitend 
    702             inrc = jstp - nit000 + 2 
    703             WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc) 
    704          END DO 
    705       ENDIF 
    706  
    707 1997  FORMAT(10X,'Time step',5X,'Sea surface temperature') 
    708 1998  FORMAT(10X,'---------',5X,'-----------------') 
    709 1999  FORMAT(10X,I9,5X,I17) 
    710        
    711    END SUBROUTINE obs_pre_sst 
    712  
    713    SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea ) 
    714       !!---------------------------------------------------------------------- 
    715       !!                    ***  ROUTINE obs_pre_seaice  *** 
    716       !! 
    717       !! ** Purpose : First level check and screening of Sea Ice observations 
    718       !! 
    719       !! ** Method  : First level check and screening of Sea Ice observations 
    720       !! 
    721       !! ** Action  :  
    722       !! 
    723       !! References : 
    724       !!    
    725       !! History : 
    726       !!        !  2007-11 (D. Lea) based on obs_pre_sst 
    727       !!---------------------------------------------------------------------- 
    728       !! * Modules used 
    729       USE domstp              ! Domain: set the time-step 
    730       USE par_oce             ! Ocean parameters 
    731       USE dom_oce, ONLY : &   ! Geographical information 
    732          & glamt,   & 
    733          & gphit,   & 
    734          & tmask,   & 
    735          & nproc 
    736       !! * Arguments 
    737       TYPE(obs_surf), INTENT(INOUT) :: seaicedata     ! Full set of Sea Ice data 
    738       TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc    ! Subset of sea ice data not failing screening 
    739       LOGICAL :: ld_seaice     ! Switch for sea ice data 
    740       LOGICAL :: ld_nea        ! Switch for rejecting observation near land 
    741       !! * Local declarations 
    742       INTEGER :: iyea0         ! Initial date 
    743       INTEGER :: imon0         !  - (year, month, day, hour, minute) 
    744       INTEGER :: iday0     
    745       INTEGER :: ihou0     
    746       INTEGER :: imin0 
    747       INTEGER :: icycle       ! Current assimilation cycle 
    748                               ! Counters for observations that 
    749       INTEGER :: iotdobs      !  - outside time domain 
    750       INTEGER :: iosdsobs     !  - outside space domain 
    751       INTEGER :: ilansobs     !  - within a model land cell 
    752       INTEGER :: inlasobs     !  - close to land 
    753       INTEGER :: igrdobs      !  - fail the grid search 
    754                               ! Global counters for observations that 
    755       INTEGER :: iotdobsmpp   !  - outside time domain 
    756       INTEGER :: iosdsobsmpp  !  - outside space domain 
    757       INTEGER :: ilansobsmpp  !  - within a model land cell 
    758       INTEGER :: inlasobsmpp  !  - close to land 
    759       INTEGER :: igrdobsmpp   !  - fail the grid search 
    760       LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    761          & llvalid            ! data selection 
    762       INTEGER :: jobs         ! Obs. loop variable 
    763       INTEGER :: jstp         ! Time loop variable 
    764       INTEGER :: inrc         ! Time index variable 
    765  
    766       IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...' 
    767  
    768       ! Initial date initialization (year, month, day, hour, minute) 
    769       iyea0 =   ndate0 / 10000 
    770       imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
    771       iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
    772       ihou0 = 0 
    773       imin0 = 0 
    774  
    775       icycle = no     ! Assimilation cycle 
    776  
    777       ! Diagnotics counters for various failures. 
    778  
    779       iotdobs  = 0 
    780       igrdobs  = 0 
    781       iosdsobs = 0 
    782       ilansobs = 0 
    783       inlasobs = 0 
    784  
    785       ! ----------------------------------------------------------------------- 
    786       ! Find time coordinate for sea ice data 
    787       ! ----------------------------------------------------------------------- 
    788  
    789       CALL obs_coo_tim( icycle, & 
    790          &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    791          &              seaicedata%nsurf,   seaicedata%nyea, seaicedata%nmon, & 
    792          &              seaicedata%nday,    seaicedata%nhou, seaicedata%nmin, & 
    793          &              seaicedata%nqc,     seaicedata%mstp, iotdobs        ) 
    794       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
    795       ! ----------------------------------------------------------------------- 
    796       ! Check for sea ice data failing the grid search 
    797       ! ----------------------------------------------------------------------- 
    798  
    799       CALL obs_coo_grd( seaicedata%nsurf,   seaicedata%mi, seaicedata%mj, & 
    800          &              seaicedata%nqc,     igrdobs                         ) 
    801       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    802  
    803       ! ----------------------------------------------------------------------- 
    804       ! Check for land points.  
    805       ! ----------------------------------------------------------------------- 
    806  
    807       CALL obs_coo_spc_2d( seaicedata%nsurf,                 & 
    808          &                 jpi,             jpj,             & 
    809          &                 seaicedata%mi,   seaicedata%mj,   &  
    810          &                 seaicedata%rlam, seaicedata%rphi, & 
    811          &                 glamt,           gphit,           & 
    812          &                 tmask(:,:,1),    seaicedata%nqc,  & 
    813          &                 iosdsobs,        ilansobs,        & 
    814          &                 inlasobs,        ld_nea           ) 
    815  
    816       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    817       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    818       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
    819  
    820       ! ----------------------------------------------------------------------- 
    821       ! Copy useful data from the seaicedata data structure to 
    822       ! the seaicedatqc data structure  
    823       ! ----------------------------------------------------------------------- 
    824  
    825       ! Allocate the selection arrays 
    826  
    827       ALLOCATE( llvalid(seaicedata%nsurf) ) 
    828        
    829       ! We want all data which has qc flags <= 0 
    830  
    831       llvalid(:)  = ( seaicedata%nqc(:)  <= 10 ) 
    832  
    833       ! The actual copying 
    834  
    835       CALL obs_surf_compress( seaicedata,     seaicedatqc,       .TRUE.,  numout, & 
    836          &                    lvalid=llvalid ) 
    837  
    838       ! Dellocate the selection arrays 
    839       DEALLOCATE( llvalid ) 
    840  
    841       ! ----------------------------------------------------------------------- 
    842       ! Print information about what observations are left after qc 
    843       ! ----------------------------------------------------------------------- 
    844  
    845       ! Update the total observation counter array 
    846        
    847       IF(lwp) THEN 
    848          WRITE(numout,*) 
    849          WRITE(numout,*) 'obs_pre_seaice :' 
    850          WRITE(numout,*) '~~~~~~~~~~~' 
    851          WRITE(numout,*) 
    852          WRITE(numout,*) ' Sea ice data outside time domain                  = ', & 
    853             &            iotdobsmpp 
    854          WRITE(numout,*) ' Remaining sea ice data that failed grid search    = ', & 
    855             &            igrdobsmpp 
    856          WRITE(numout,*) ' Remaining sea ice data outside space domain       = ', & 
    857             &            iosdsobsmpp 
    858          WRITE(numout,*) ' Remaining sea ice data at land points             = ', & 
    859             &            ilansobsmpp 
    860          IF (ld_nea) THEN 
    861             WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', & 
    862                &            inlasobsmpp 
    863          ELSE 
    864             WRITE(numout,*) ' Remaining sea ice data near land points (kept)    = ', & 
    865                &            inlasobsmpp 
    866          ENDIF 
    867          WRITE(numout,*) ' Sea ice data accepted                             = ', & 
    868             &            seaicedatqc%nsurfmpp 
    869  
    870          WRITE(numout,*) 
    871          WRITE(numout,*) ' Number of observations per time step :' 
    872          WRITE(numout,*) 
    873          WRITE(numout,1997) 
    874          WRITE(numout,1998) 
    875       ENDIF 
    876        
    877       DO jobs = 1, seaicedatqc%nsurf 
    878          inrc = seaicedatqc%mstp(jobs) + 2 - nit000 
    879          seaicedatqc%nsstp(inrc)  = seaicedatqc%nsstp(inrc) + 1 
    880       END DO 
    881        
    882       CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, & 
    883          &                       nitend - nit000 + 2 ) 
    884  
    885       IF ( lwp ) THEN 
    886          DO jstp = nit000 - 1, nitend 
    887             inrc = jstp - nit000 + 2 
    888             WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc) 
    889          END DO 
    890       ENDIF 
    891  
    892 1997  FORMAT(10X,'Time step',5X,'Sea ice data           ') 
    893 1998  FORMAT(10X,'---------',5X,'-----------------') 
    894 1999  FORMAT(10X,I9,5X,I17) 
    895        
    896    END SUBROUTINE obs_pre_seaice 
    897  
    898    SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav ) 
    899       !!---------------------------------------------------------------------- 
    900       !!                    ***  ROUTINE obs_pre_taovel  *** 
    901       !! 
    902       !! ** Purpose : First level check and screening of U and V profiles 
    903       !! 
    904       !! ** Method  : First level check and screening of U and V profiles 
     234   END SUBROUTINE obs_pre_surf 
     235 
     236 
     237   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav ) 
     238      !!---------------------------------------------------------------------- 
     239      !!                    ***  ROUTINE obs_pre_prof  *** 
     240      !! 
     241      !! ** Purpose : First level check and screening of profiles 
     242      !! 
     243      !! ** Method  : First level check and screening of profiles 
    905244      !! 
    906245      !! History : 
     
    908247      !!        !  2008-09  (M. Valdivieso) : TAO velocity data 
    909248      !!        !  2009-01  (K. Mogensen) : New feedback strictuer 
     249      !!        !  2015-02  (M. Martin) : Combined profile routine. 
    910250      !! 
    911251      !!---------------------------------------------------------------------- 
     
    962302      INTEGER :: inrc         ! Time index variable 
    963303 
    964       IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data' 
     304      IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data' 
    965305 
    966306      ! Initial date initialization (year, month, day, hour, minute) 
     
    1186526999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 
    1187527 
    1188    END SUBROUTINE obs_pre_vel 
     528   END SUBROUTINE obs_pre_prof 
    1189529 
    1190530   SUBROUTINE obs_coo_tim( kcycle, & 
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r4990 r5659  
    3333   PRIVATE 
    3434 
    35    PUBLIC obs_rea_pro_dri  ! Read the profile observations  
     35   PUBLIC obs_rea_prof  ! Read the profile observations  
    3636 
    3737   !!---------------------------------------------------------------------- 
     
    4343CONTAINS 
    4444  
    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 ) 
     45   SUBROUTINE obs_rea_prof( profdata, knumfiles, cfilenames, & 
     46      &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
     47      &                     ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, & 
     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  :  
     
    6766    
    6867      !! * Arguments 
    69       INTEGER ::  kformat    ! Format of input data 
    70       !                      ! 1: ENACT 
    71       !                      ! 2: Coriolis 
    7268      TYPE(obs_prof), INTENT(OUT) ::  profdata     ! Profile data to be read 
    7369      INTEGER, INTENT(IN) :: knumfiles      ! Number of files to read in 
     
    8985 
    9086      !! * Local declarations 
    91       CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_pro_dri' 
     87      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    9288      INTEGER :: jvar 
    9389      INTEGER :: ji 
     
    195191         !--------------------------------------------------------------------- 
    196192          
    197          iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, & 
     193         iflag = nf90_open( TRIM( cfilenames(jj) ), nf90_nowrite, & 
    198194            &                      i_file_id ) 
    199195          
     
    202198            IF ( ldignmis ) THEN 
    203199               inpfiles(jj)%nobs = 0 
    204                CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // & 
     200               CALL ctl_warn( 'File ' // TRIM( cfilenames(jj) ) // & 
    205201                  &           ' not found' ) 
    206202            ELSE  
    207                CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // & 
     203               CALL ctl_stop( 'File ' // TRIM( cfilenames(jj) ) // & 
    208204                  &           ' not found' ) 
    209205            ENDIF 
     
    220216            !  Read the profile file into inpfiles 
    221217            !------------------------------------------------------------------ 
    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. ) 
    253             ELSE 
    254                CALL ctl_stop( 'File format unknown' ) 
     218            CALL init_obfbdata( inpfiles(jj) ) 
     219            IF(lwp) THEN 
     220               WRITE(numout,*) 
     221               WRITE(numout,*)'Reading from feedback file :', & 
     222                  &           TRIM( cfilenames(jj) ) 
     223            ENDIF 
     224            CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), & 
     225               &                ldgrid = .TRUE. ) 
     226                
     227            IF ( inpfiles(jj)%nvar < 2 ) THEN 
     228               CALL ctl_stop( 'Feedback format error' ) 
     229               RETURN 
     230            ENDIF 
     231            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 
     232               CALL ctl_stop( 'Model not in input data' ) 
     233               RETURN 
    255234            ENDIF 
    256235             
     
    836815      DEALLOCATE( inpfiles ) 
    837816 
    838    END SUBROUTINE obs_rea_pro_dri 
     817   END SUBROUTINE obs_rea_prof 
    839818 
    840819END MODULE obs_read_prof 
  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r4990 r5659  
    66 
    77   !!---------------------------------------------------------------------- 
    8    !!   obs_wri_p3d   : Write profile observation diagnostics in NetCDF format 
    9    !!   obs_wri_sla   : Write SLA observation related diagnostics 
    10    !!   obs_wri_sst   : Write SST observation related diagnostics 
    11    !!   obs_wri_seaice: Write seaice observation related diagnostics 
    12    !!   obs_wri_vel   : Write velocity observation diagnostics in NetCDF format 
     8   !!   obs_wri_prof   : Write profile observations in feedback format 
     9   !!   obs_wri_surf   : Write surface observations in feedback format 
    1310   !!   obs_wri_stats : Print basic statistics on the data being written out 
    1411   !!---------------------------------------------------------------------- 
     
    3027   USE obs_conv             ! Conversion between units 
    3128   USE obs_const 
    32    USE obs_sla_types 
    3329   USE obs_rot_vel          ! Rotation of velocities 
    3430   USE obs_mpp              ! MPP support routines for observation diagnostics 
     
    3935   !! * Routine accessibility 
    4036   PRIVATE 
    41    PUBLIC obs_wri_p3d, &    ! Write profile observation related diagnostics 
    42       &   obs_wri_sla, &    ! Write SLA observation related diagnostics 
    43       &   obs_wri_sst, &    ! Write SST observation related diagnostics 
    44       &   obs_wri_sss, &    ! Write SSS observation related diagnostics 
    45       &   obs_wri_seaice, & ! Write seaice observation related diagnostics 
    46       &   obs_wri_vel, &    ! Write velocity observation related diagnostics 
     37   PUBLIC obs_wri_prof, &    ! Write profile observation files 
     38      &   obs_wri_surf, &    ! Write surface observation files 
    4739      &   obswriinfo 
    4840    
     
    6355CONTAINS 
    6456 
    65    SUBROUTINE obs_wri_p3d( cprefix, profdata, padd, pext ) 
     57   SUBROUTINE obs_wri_prof( cobstype, profdata, padd, pext ) 
    6658      !!----------------------------------------------------------------------- 
    6759      !! 
    68       !!                     *** ROUTINE obs_wri_p3d  *** 
    69       !! 
    70       !! ** Purpose : Write temperature and salinity (profile) observation  
    71       !!              related diagnostics 
     60      !!                     *** ROUTINE obs_wri_prof  *** 
     61      !! 
     62      !! ** Purpose : Write profile feedback files 
    7263      !! 
    7364      !! ** Method  : NetCDF 
     
    8273      !!      ! 07-03  (K. Mogensen) General handling of profiles 
    8374      !!      ! 09-01  (K. Mogensen) New feedback format 
     75      !!      ! 15-02  (M. Martin) Combined routine for writing profiles 
    8476      !!----------------------------------------------------------------------- 
    8577 
     
    8779 
    8880      !! * Arguments 
    89       CHARACTER(LEN=*), INTENT(IN) :: cprefix        ! Prefix for output files 
     81      CHARACTER(LEN=*), INTENT(IN) :: cobstype        ! Prefix for output files 
    9082      TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data 
    9183      TYPE(obswriinfo), OPTIONAL :: padd             ! Additional info for each variable 
     
    125117         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 
    126118      END DO 
    127       CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 
    128          &                 1 + nadd, 1 + next, .TRUE. ) 
    129  
    130       fbdata%cname(1)      = 'POTM' 
    131       fbdata%cname(2)      = 'PSAL' 
    132       fbdata%coblong(1)    = 'Potential temperature' 
    133       fbdata%coblong(2)    = 'Practical salinity' 
    134       fbdata%cobunit(1)    = 'Degrees centigrade' 
    135       fbdata%cobunit(2)    = 'PSU' 
    136       fbdata%cextname(1)   = 'TEMP' 
    137       fbdata%cextlong(1)   = 'Insitu temperature' 
    138       fbdata%cextunit(1)   = 'Degrees centigrade' 
    139       DO je = 1, next 
    140          fbdata%cextname(1+je) = pext%cdname(je) 
    141          fbdata%cextlong(1+je) = pext%cdlong(je,1) 
    142          fbdata%cextunit(1+je) = pext%cdunit(je,1) 
    143       END DO 
     119 
     120      SELECT CASE ( TRIM(cobstype) ) 
     121      CASE('prof') 
     122 
     123         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 
     124            &                 1 + nadd, 1 + next, .TRUE. ) 
     125         fbdata%cname(1)      = 'POTM' 
     126         fbdata%cname(2)      = 'PSAL' 
     127         fbdata%coblong(1)    = 'Potential temperature' 
     128         fbdata%coblong(2)    = 'Practical salinity' 
     129         fbdata%cobunit(1)    = 'Degrees centigrade' 
     130         fbdata%cobunit(2)    = 'PSU' 
     131         fbdata%cextname(1)   = 'TEMP' 
     132         fbdata%cextlong(1)   = 'Insitu temperature' 
     133         fbdata%cextunit(1)   = 'Degrees centigrade' 
     134         fbdata%caddlong(1,1) = 'Model interpolated potential temperature' 
     135         fbdata%caddlong(1,2) = 'Model interpolated practical salinity' 
     136         fbdata%caddunit(1,1) = 'Degrees centigrade' 
     137         fbdata%caddunit(1,2) = 'PSU' 
     138         fbdata%cgrid(:)      = 'T' 
     139         DO je = 1, next 
     140            fbdata%cextname(1+je) = pext%cdname(je) 
     141            fbdata%cextlong(1+je) = pext%cdlong(je,1) 
     142            fbdata%cextunit(1+je) = pext%cdunit(je,1) 
     143         END DO 
     144         DO ja = 1, nadd 
     145            fbdata%caddname(1+ja) = padd%cdname(ja) 
     146            DO jvar = 1, 2 
     147               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar) 
     148               fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar) 
     149            END DO 
     150         END DO 
     151       
     152      CASE('vel') 
     153       
     154         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. ) 
     155         fbdata%cname(1)      = 'UVEL' 
     156         fbdata%cname(2)      = 'VVEL' 
     157         fbdata%coblong(1)    = 'Zonal velocity' 
     158         fbdata%coblong(2)    = 'Meridional velocity' 
     159         fbdata%cobunit(1)    = 'm/s' 
     160         fbdata%cobunit(2)    = 'm/s' 
     161         DO je = 1, next 
     162            fbdata%cextname(je) = pext%cdname(je) 
     163            fbdata%cextlong(je) = pext%cdlong(je,1) 
     164            fbdata%cextunit(je) = pext%cdunit(je,1) 
     165         END DO 
     166         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity' 
     167         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity' 
     168         fbdata%caddunit(1,1) = 'm/s' 
     169         fbdata%caddunit(1,2) = 'm/s' 
     170         fbdata%caddname(2)   = 'HxG' 
     171         fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)' 
     172         fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)' 
     173         fbdata%caddunit(2,1) = 'm/s' 
     174         fbdata%caddunit(2,2) = 'm/s'  
     175         fbdata%cgrid(1)      = 'U'  
     176         fbdata%cgrid(2)      = 'V' 
     177         DO ja = 1, nadd 
     178            fbdata%caddname(2+ja) = padd%cdname(ja) 
     179            fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 
     180            fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 
     181         END DO 
     182         ALLOCATE( & 
     183            & zu(profdata%nvprot(1)), & 
     184            & zv(profdata%nvprot(2))  & 
     185            & ) 
     186         CALL obs_rotvel( profdata, k2dint, zu, zv ) 
     187 
     188      END SELECT 
     189          
    144190      fbdata%caddname(1)   = 'Hx' 
    145       fbdata%caddlong(1,1) = 'Model interpolated potential temperature' 
    146       fbdata%caddlong(1,2) = 'Model interpolated practical salinity' 
    147       fbdata%caddunit(1,1) = 'Degrees centigrade' 
    148       fbdata%caddunit(1,2) = 'PSU' 
    149       fbdata%cgrid(:)      = 'T' 
    150       DO ja = 1, nadd 
    151          fbdata%caddname(1+ja) = padd%cdname(ja) 
    152          DO jvar = 1, 2 
    153             fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar) 
    154             fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar) 
    155          END DO 
    156       END DO 
    157191          
    158       WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     192      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc 
    159193 
    160194      IF(lwp) THEN 
    161195         WRITE(numout,*) 
    162          WRITE(numout,*)'obs_wri_p3d :' 
     196         WRITE(numout,*)'obs_wri_prof :' 
    163197         WRITE(numout,*)'~~~~~~~~~~~~~' 
    164198         WRITE(numout,*)'Writing profile feedback file : ',TRIM(cfname) 
     
    208242            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    209243               ik = profdata%var(jvar)%nvlidx(jk) 
    210                fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
     244               IF ( TRIM(cobstype) == 'prof' ) THEN 
     245                  fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
     246               ELSE IF ( TRIM(cobstype) == 'vel' ) THEN 
     247                  IF ( jvar == 1 ) THEN 
     248                     fbdata%padd(ik,jo,1,jvar) = zu(jk) 
     249                  ELSE 
     250                     fbdata%padd(ik,jo,1,jvar) = zv(jk) 
     251                  ENDIF 
     252                  fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk) 
     253               ENDIF 
    211254               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk) 
    212255               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk) 
     
    237280      END DO 
    238281 
    239       ! Convert insitu temperature to potential temperature using the model 
    240       ! salinity if no potential temperature 
    241       DO jo = 1, fbdata%nobs 
    242          IF ( fbdata%pphi(jo) < 9999.0 ) THEN 
    243             DO jk = 1, fbdata%nlev 
    244                IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. & 
    245                   & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. & 
    246                   & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. & 
    247                   & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN 
    248                   zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), & 
    249                      &              REAL(fbdata%pphi(jo),wp) ) 
    250                   fbdata%pob(jk,jo,1) = potemp( & 
    251                      &                     REAL(fbdata%padd(jk,jo,1,2), wp),  & 
    252                      &                     REAL(fbdata%pext(jk,jo,1), wp), & 
    253                      &                     zpres, 0.0_wp ) 
    254                ENDIF 
    255             END DO 
    256          ENDIF 
    257       END DO 
     282      IF ( TRIM(cobstype) == 'prof' ) THEN 
     283         ! Convert insitu temperature to potential temperature using the model 
     284         ! salinity if no potential temperature 
     285         DO jo = 1, fbdata%nobs 
     286            IF ( fbdata%pphi(jo) < 9999.0 ) THEN 
     287               DO jk = 1, fbdata%nlev 
     288                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. & 
     289                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. & 
     290                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. & 
     291                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN 
     292                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), & 
     293                        &              REAL(fbdata%pphi(jo),wp) ) 
     294                     fbdata%pob(jk,jo,1) = potemp( & 
     295                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  & 
     296                        &                     REAL(fbdata%pext(jk,jo,1), wp), & 
     297                        &                     zpres, 0.0_wp ) 
     298                  ENDIF 
     299               END DO 
     300            ENDIF 
     301         END DO 
     302      ENDIF 
    258303       
    259304      ! Write the obfbdata structure 
     
    265310      CALL dealloc_obfbdata( fbdata ) 
    266311      
    267    END SUBROUTINE obs_wri_p3d 
    268  
    269    SUBROUTINE obs_wri_sla( cprefix, sladata, padd, pext ) 
     312   END SUBROUTINE obs_wri_prof 
     313 
     314   SUBROUTINE obs_wri_surf( cobstype, surfdata, padd, pext ) 
    270315      !!----------------------------------------------------------------------- 
    271316      !! 
    272       !!                     *** ROUTINE obs_wri_sla  *** 
    273       !! 
    274       !! ** Purpose : Write SLA observation diagnostics 
    275       !!              related  
     317      !!                     *** ROUTINE obs_wri_surf  *** 
     318      !! 
     319      !! ** Purpose : Write surface observation files 
    276320      !! 
    277321      !! ** Method  : NetCDF 
     
    281325      !!      ! 07-03  (K. Mogensen) Original 
    282326      !!      ! 09-01  (K. Mogensen) New feedback format. 
     327      !!      ! 15-02  (M. Martin) Combined surface writing routine. 
    283328      !!----------------------------------------------------------------------- 
    284329 
     
    287332 
    288333      !! * Arguments 
    289       CHARACTER(LEN=*), INTENT(IN) :: cprefix          ! Prefix for output files 
    290       TYPE(obs_surf), INTENT(INOUT) :: sladata         ! Full set of SLAa 
     334      CHARACTER(LEN=*), INTENT(IN) :: cobstype          ! Prefix for output files 
     335      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data 
    291336      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable 
    292337      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info 
     
    295340      TYPE(obfbdata) :: fbdata 
    296341      CHARACTER(LEN=40) :: cfname         ! netCDF filename 
    297       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_sla' 
     342      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    298343      INTEGER :: jo 
    299344      INTEGER :: ja 
     
    316361      CALL init_obfbdata( fbdata ) 
    317362 
    318       CALL alloc_obfbdata( fbdata, 1, sladata%nsurf, 1, & 
     363      CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    319364         &                 2 + nadd, 1 + next, .TRUE. ) 
    320365 
    321       fbdata%cname(1)      = 'SLA' 
    322       fbdata%coblong(1)    = 'Sea level anomaly' 
    323       fbdata%cobunit(1)    = 'Metres' 
    324       fbdata%cextname(1)   = 'MDT' 
    325       fbdata%cextlong(1)   = 'Mean dynamic topography' 
    326       fbdata%cextunit(1)   = 'Metres' 
    327       DO je = 1, next 
    328          fbdata%cextname(1+je) = pext%cdname(je) 
    329          fbdata%cextlong(1+je) = pext%cdlong(je,1) 
    330          fbdata%cextunit(1+je) = pext%cdunit(je,1) 
    331       END DO 
     366      SELECT CASE ( TRIM(cobstype) ) 
     367      CASE('sla') 
     368 
     369         fbdata%cname(1)      = 'SLA' 
     370         fbdata%coblong(1)    = 'Sea level anomaly' 
     371         fbdata%cobunit(1)    = 'Metres' 
     372         fbdata%cextname(1)   = 'MDT' 
     373         fbdata%cextlong(1)   = 'Mean dynamic topography' 
     374         fbdata%cextunit(1)   = 'Metres' 
     375         DO je = 1, next 
     376            fbdata%cextname(je) = pext%cdname(je) 
     377            fbdata%cextlong(je) = pext%cdlong(je,1) 
     378            fbdata%cextunit(je) = pext%cdunit(je,1) 
     379         END DO 
     380         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT' 
     381         fbdata%caddunit(1,1) = 'Metres'  
     382         fbdata%caddname(2)   = 'SSH' 
     383         fbdata%caddlong(2,1) = 'Model Sea surface height' 
     384         fbdata%caddunit(2,1) = 'Metres' 
     385         fbdata%cgrid(1)      = 'T' 
     386         DO ja = 1, nadd 
     387            fbdata%caddname(2+ja) = padd%cdname(ja) 
     388            fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 
     389            fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 
     390         END DO 
     391 
     392      CASE('sst') 
     393 
     394         fbdata%cname(1)      = 'SST' 
     395         fbdata%coblong(1)    = 'Sea surface temperature' 
     396         fbdata%cobunit(1)    = 'Degree centigrade' 
     397         DO je = 1, next 
     398            fbdata%cextname(je) = pext%cdname(je) 
     399            fbdata%cextlong(je) = pext%cdlong(je,1) 
     400            fbdata%cextunit(je) = pext%cdunit(je,1) 
     401         END DO 
     402         fbdata%caddlong(1,1) = 'Model interpolated SST' 
     403         fbdata%caddunit(1,1) = 'Degree centigrade' 
     404         fbdata%cgrid(1)      = 'T' 
     405         DO ja = 1, nadd 
     406            fbdata%caddname(1+ja) = padd%cdname(ja) 
     407            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     408            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     409         END DO 
     410 
     411      CASE('sss') 
     412 
     413         fbdata%cname(1)      = 'SSS' 
     414         fbdata%coblong(1)    = 'Sea surface salinity' 
     415         fbdata%cobunit(1)    = 'PSU' 
     416         DO je = 1, next 
     417            fbdata%cextname(je) = pext%cdname(je) 
     418            fbdata%cextlong(je) = pext%cdlong(je,1) 
     419            fbdata%cextunit(je) = pext%cdunit(je,1) 
     420         END DO 
     421         fbdata%caddlong(1,1) = 'Model interpolated SSS' 
     422         fbdata%caddunit(1,1) = 'PSU' 
     423         fbdata%cgrid(1)      = 'T' 
     424         DO ja = 1, nadd 
     425            fbdata%caddname(1+ja) = padd%cdname(ja) 
     426            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     427            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     428         END DO 
     429 
     430      CASE('seaice') 
     431 
     432         fbdata%cname(1)      = 'SEAICE' 
     433         fbdata%coblong(1)    = 'Sea ice' 
     434         fbdata%cobunit(1)    = 'Fraction' 
     435         DO je = 1, next 
     436            fbdata%cextname(je) = pext%cdname(je) 
     437            fbdata%cextlong(je) = pext%cdlong(je,1) 
     438            fbdata%cextunit(je) = pext%cdunit(je,1) 
     439         END DO 
     440         fbdata%caddlong(1,1) = 'Model interpolated ICE' 
     441         fbdata%caddunit(1,1) = 'Fraction' 
     442         fbdata%cgrid(1)      = 'T' 
     443         DO ja = 1, nadd 
     444            fbdata%caddname(1+ja) = padd%cdname(ja) 
     445            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     446            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     447         END DO 
     448 
     449      END SELECT 
     450 
    332451      fbdata%caddname(1)   = 'Hx' 
    333       fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT' 
    334       fbdata%caddunit(1,1) = 'Metres'  
    335       fbdata%caddname(2)   = 'SSH' 
    336       fbdata%caddlong(2,1) = 'Model Sea surface height' 
    337       fbdata%caddunit(2,1) = 'Metres' 
    338       fbdata%cgrid(1)      = 'T' 
    339       DO ja = 1, nadd 
    340          fbdata%caddname(2+ja) = padd%cdname(ja) 
    341          fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 
    342          fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 
    343       END DO 
    344  
    345       WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     452 
     453      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc 
    346454 
    347455      IF(lwp) THEN 
    348456         WRITE(numout,*) 
    349          WRITE(numout,*)'obs_wri_sla :' 
     457         WRITE(numout,*)'obs_wri_surf :' 
    350458         WRITE(numout,*)'~~~~~~~~~~~~~' 
    351          WRITE(numout,*)'Writing SLA feedback file : ',TRIM(cfname) 
     459         WRITE(numout,*)'Writing surface feedback file : ',TRIM(cfname) 
    352460      ENDIF 
    353461 
    354462      ! Transform obs_prof data structure into obfbdata structure 
    355463      fbdata%cdjuldref = '19500101000000' 
    356       DO jo = 1, sladata%nsurf 
    357          fbdata%plam(jo)      = sladata%rlam(jo) 
    358          fbdata%pphi(jo)      = sladata%rphi(jo) 
    359          WRITE(fbdata%cdtyp(jo),'(I4)') sladata%ntyp(jo) 
     464      DO jo = 1, surfdata%nsurf 
     465         fbdata%plam(jo)      = surfdata%rlam(jo) 
     466         fbdata%pphi(jo)      = surfdata%rphi(jo) 
     467         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo) 
    360468         fbdata%ivqc(jo,:)    = 0 
    361469         fbdata%ivqcf(:,jo,:) = 0 
    362          IF ( sladata%nqc(jo) > 10 ) THEN 
     470         IF ( surfdata%nqc(jo) > 10 ) THEN 
    363471            fbdata%ioqc(jo)    = 4 
    364472            fbdata%ioqcf(1,jo) = 0 
    365             fbdata%ioqcf(2,jo) = sladata%nqc(jo) - 10 
     473            fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 
    366474         ELSE 
    367             fbdata%ioqc(jo)    = sladata%nqc(jo) 
     475            fbdata%ioqc(jo)    = surfdata%nqc(jo) 
    368476            fbdata%ioqcf(:,jo) = 0 
    369477         ENDIF 
     
    372480         fbdata%itqc(jo)      = 0 
    373481         fbdata%itqcf(:,jo)   = 0 
    374          fbdata%cdwmo(jo)     = sladata%cwmo(jo) 
    375          fbdata%kindex(jo)    = sladata%nsfil(jo) 
     482         fbdata%cdwmo(jo)     = surfdata%cwmo(jo) 
     483         fbdata%kindex(jo)    = surfdata%nsfil(jo) 
    376484         IF (ln_grid_global) THEN 
    377             fbdata%iobsi(jo,1) = sladata%mi(jo) 
    378             fbdata%iobsj(jo,1) = sladata%mj(jo) 
     485            fbdata%iobsi(jo,1) = surfdata%mi(jo) 
     486            fbdata%iobsj(jo,1) = surfdata%mj(jo) 
    379487         ELSE 
    380             fbdata%iobsi(jo,1) = mig(sladata%mi(jo)) 
    381             fbdata%iobsj(jo,1) = mjg(sladata%mj(jo)) 
     488            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo)) 
     489            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo)) 
    382490         ENDIF 
    383491         CALL greg2jul( 0, & 
    384             &           sladata%nmin(jo), & 
    385             &           sladata%nhou(jo), & 
    386             &           sladata%nday(jo), & 
    387             &           sladata%nmon(jo), & 
    388             &           sladata%nyea(jo), & 
     492            &           surfdata%nmin(jo), & 
     493            &           surfdata%nhou(jo), & 
     494            &           surfdata%nday(jo), & 
     495            &           surfdata%nmon(jo), & 
     496            &           surfdata%nyea(jo), & 
    389497            &           fbdata%ptim(jo),   & 
    390498            &           krefdate = 19500101 ) 
    391          fbdata%padd(1,jo,1,1) = sladata%rmod(jo,1) 
    392          fbdata%padd(1,jo,2,1) = sladata%rext(jo,1) 
    393          fbdata%pob(1,jo,1)    = sladata%robs(jo,1)  
     499         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1) 
     500         IF ( TRIM(cobstype) == 'sla' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1) 
     501         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1)  
    394502         fbdata%pdep(1,jo)     = 0.0 
    395503         fbdata%idqc(1,jo)     = 0 
    396504         fbdata%idqcf(:,1,jo)  = 0 
    397          IF ( sladata%nqc(jo) > 10 ) THEN 
     505         IF ( surfdata%nqc(jo) > 10 ) THEN 
    398506            fbdata%ivqc(jo,1)       = 4 
    399507            fbdata%ivlqc(1,jo,1)    = 4 
    400508            fbdata%ivlqcf(1,1,jo,1) = 0 
    401             fbdata%ivlqcf(2,1,jo,1) = sladata%nqc(jo) - 10 
     509            fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 
    402510         ELSE 
    403             fbdata%ivqc(jo,1)       = sladata%nqc(jo) 
    404             fbdata%ivlqc(1,jo,1)    = sladata%nqc(jo) 
     511            fbdata%ivqc(jo,1)       = surfdata%nqc(jo) 
     512            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo) 
    405513            fbdata%ivlqcf(:,1,jo,1) = 0 
    406514         ENDIF 
    407515         fbdata%iobsk(1,jo,1)  = 0 
    408          fbdata%pext(1,jo,1) = sladata%rext(jo,2) 
     516         IF ( TRIM(cobstype) == 'sla' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 
    409517         DO ja = 1, nadd 
    410518            fbdata%padd(1,jo,2+ja,1) = & 
    411                & sladata%rext(jo,padd%ipoint(ja)) 
     519               & surfdata%rext(jo,padd%ipoint(ja)) 
    412520         END DO 
    413521         DO je = 1, next 
    414522            fbdata%pext(1,jo,1+je) = & 
    415                & sladata%rext(jo,pext%ipoint(je)) 
     523               & surfdata%rext(jo,pext%ipoint(je)) 
    416524         END DO 
    417525      END DO 
     
    425533      CALL dealloc_obfbdata( fbdata ) 
    426534 
    427    END SUBROUTINE obs_wri_sla 
    428  
    429    SUBROUTINE obs_wri_sst( cprefix, sstdata, padd, pext ) 
    430       !!----------------------------------------------------------------------- 
    431       !! 
    432       !!                     *** ROUTINE obs_wri_sst  *** 
    433       !! 
    434       !! ** Purpose : Write SST observation diagnostics 
    435       !!              related  
    436       !! 
    437       !! ** Method  : NetCDF 
    438       !!  
    439       !! ** Action  : 
    440       !! 
    441       !!      ! 07-07  (S. Ricci) Original 
    442       !!      ! 09-01  (K. Mogensen) New feedback format. 
    443       !!----------------------------------------------------------------------- 
    444  
    445       !! * Modules used 
    446       IMPLICIT NONE 
    447  
    448       !! * Arguments 
    449       CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
    450       TYPE(obs_surf), INTENT(INOUT) :: sstdata      ! Full set of SST 
    451       TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
    452       TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
    453  
    454       !! * Local declarations  
    455       TYPE(obfbdata) :: fbdata 
    456       CHARACTER(LEN=40) ::  cfname             ! netCDF filename 
    457       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_sst' 
    458       INTEGER :: jo 
    459       INTEGER :: ja 
    460       INTEGER :: je 
    461       INTEGER :: nadd 
    462       INTEGER :: next 
    463  
    464       IF ( PRESENT( padd ) ) THEN 
    465          nadd = padd%inum 
    466       ELSE 
    467          nadd = 0 
    468       ENDIF 
    469  
    470       IF ( PRESENT( pext ) ) THEN 
    471          next = pext%inum 
    472       ELSE 
    473          next = 0 
    474       ENDIF 
    475  
    476       CALL init_obfbdata( fbdata ) 
    477  
    478       CALL alloc_obfbdata( fbdata, 1, sstdata%nsurf, 1, & 
    479          &                 1 + nadd, next, .TRUE. ) 
    480  
    481       fbdata%cname(1)      = 'SST' 
    482       fbdata%coblong(1)    = 'Sea surface temperature' 
    483       fbdata%cobunit(1)    = 'Degree centigrade' 
    484       DO je = 1, next 
    485          fbdata%cextname(je) = pext%cdname(je) 
    486          fbdata%cextlong(je) = pext%cdlong(je,1) 
    487          fbdata%cextunit(je) = pext%cdunit(je,1) 
    488       END DO 
    489       fbdata%caddname(1)   = 'Hx' 
    490       fbdata%caddlong(1,1) = 'Model interpolated SST' 
    491       fbdata%caddunit(1,1) = 'Degree centigrade' 
    492       fbdata%cgrid(1)      = 'T' 
    493       DO ja = 1, nadd 
    494          fbdata%caddname(1+ja) = padd%cdname(ja) 
    495          fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    496          fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    497       END DO 
    498  
    499       WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
    500  
    501       IF(lwp) THEN 
    502          WRITE(numout,*) 
    503          WRITE(numout,*)'obs_wri_sst :' 
    504          WRITE(numout,*)'~~~~~~~~~~~~~' 
    505          WRITE(numout,*)'Writing SST feedback file : ',TRIM(cfname) 
    506       ENDIF 
    507  
    508       ! Transform obs_prof data structure into obfbdata structure 
    509       fbdata%cdjuldref = '19500101000000' 
    510       DO jo = 1, sstdata%nsurf 
    511          fbdata%plam(jo)      = sstdata%rlam(jo) 
    512          fbdata%pphi(jo)      = sstdata%rphi(jo) 
    513          WRITE(fbdata%cdtyp(jo),'(I4)') sstdata%ntyp(jo) 
    514          fbdata%ivqc(jo,:)    = 0 
    515          fbdata%ivqcf(:,jo,:) = 0 
    516          IF ( sstdata%nqc(jo) > 10 ) THEN 
    517             fbdata%ioqc(jo)    = 4 
    518             fbdata%ioqcf(1,jo) = 0 
    519             fbdata%ioqcf(2,jo) = sstdata%nqc(jo) - 10 
    520          ELSE 
    521             fbdata%ioqc(jo)    = MAX(sstdata%nqc(jo),1) 
    522             fbdata%ioqcf(:,jo) = 0 
    523          ENDIF 
    524          fbdata%ipqc(jo)      = 0 
    525          fbdata%ipqcf(:,jo)   = 0 
    526          fbdata%itqc(jo)      = 0 
    527          fbdata%itqcf(:,jo)   = 0 
    528          fbdata%cdwmo(jo)     = '' 
    529          fbdata%kindex(jo)    = sstdata%nsfil(jo) 
    530          IF (ln_grid_global) THEN 
    531             fbdata%iobsi(jo,1) = sstdata%mi(jo) 
    532             fbdata%iobsj(jo,1) = sstdata%mj(jo) 
    533          ELSE 
    534             fbdata%iobsi(jo,1) = mig(sstdata%mi(jo)) 
    535             fbdata%iobsj(jo,1) = mjg(sstdata%mj(jo)) 
    536          ENDIF 
    537          CALL greg2jul( 0, & 
    538             &           sstdata%nmin(jo), & 
    539             &           sstdata%nhou(jo), & 
    540             &           sstdata%nday(jo), & 
    541             &           sstdata%nmon(jo), & 
    542             &           sstdata%nyea(jo), & 
    543             &           fbdata%ptim(jo),   & 
    544             &           krefdate = 19500101 ) 
    545          fbdata%padd(1,jo,1,1) = sstdata%rmod(jo,1) 
    546          fbdata%pob(1,jo,1)    = sstdata%robs(jo,1) 
    547          fbdata%pdep(1,jo)     = 0.0 
    548          fbdata%idqc(1,jo)     = 0 
    549          fbdata%idqcf(:,1,jo)  = 0 
    550          IF ( sstdata%nqc(jo) > 10 ) THEN 
    551             fbdata%ivqc(jo,1)       = 4 
    552             fbdata%ivlqc(1,jo,1)    = 4 
    553             fbdata%ivlqcf(1,1,jo,1) = 0 
    554             fbdata%ivlqcf(2,1,jo,1) = sstdata%nqc(jo) - 10 
    555          ELSE 
    556             fbdata%ivqc(jo,1)       = MAX(sstdata%nqc(jo),1) 
    557             fbdata%ivlqc(1,jo,1)    = MAX(sstdata%nqc(jo),1) 
    558             fbdata%ivlqcf(:,1,jo,1) = 0 
    559          ENDIF 
    560          fbdata%iobsk(1,jo,1)  = 0 
    561          DO ja = 1, nadd 
    562             fbdata%padd(1,jo,1+ja,1) = & 
    563                & sstdata%rext(jo,padd%ipoint(ja)) 
    564          END DO 
    565          DO je = 1, next 
    566             fbdata%pext(1,jo,je) = & 
    567                & sstdata%rext(jo,pext%ipoint(je)) 
    568          END DO 
    569  
    570       END DO 
    571  
    572       ! Write the obfbdata structure 
    573  
    574       CALL write_obfbdata( cfname, fbdata ) 
    575  
    576       ! Output some basic statistics 
    577       CALL obs_wri_stats( fbdata ) 
    578  
    579       CALL dealloc_obfbdata( fbdata ) 
    580  
    581    END SUBROUTINE obs_wri_sst 
    582  
    583    SUBROUTINE obs_wri_sss 
    584    END SUBROUTINE obs_wri_sss 
    585  
    586    SUBROUTINE obs_wri_seaice( cprefix, seaicedata, padd, pext ) 
    587       !!----------------------------------------------------------------------- 
    588       !! 
    589       !!                     *** ROUTINE obs_wri_seaice  *** 
    590       !! 
    591       !! ** Purpose : Write sea ice observation diagnostics 
    592       !!              related  
    593       !! 
    594       !! ** Method  : NetCDF 
    595       !!  
    596       !! ** Action  : 
    597       !! 
    598       !!      ! 07-07  (S. Ricci) Original 
    599       !!      ! 09-01  (K. Mogensen) New feedback format. 
    600       !!----------------------------------------------------------------------- 
    601  
    602       !! * Modules used 
    603       IMPLICIT NONE 
    604  
    605       !! * Arguments 
    606       CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
    607       TYPE(obs_surf), INTENT(INOUT) :: seaicedata   ! Full set of sea ice 
    608       TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
    609       TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
    610  
    611       !! * Local declarations  
    612       TYPE(obfbdata) :: fbdata 
    613       CHARACTER(LEN=40) :: cfname             ! netCDF filename 
    614       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_seaice' 
    615       INTEGER :: jo 
    616       INTEGER :: ja 
    617       INTEGER :: je 
    618       INTEGER :: nadd 
    619       INTEGER :: next 
    620  
    621       IF ( PRESENT( padd ) ) THEN 
    622          nadd = padd%inum 
    623       ELSE 
    624          nadd = 0 
    625       ENDIF 
    626  
    627       IF ( PRESENT( pext ) ) THEN 
    628          next = pext%inum 
    629       ELSE 
    630          next = 0 
    631       ENDIF 
    632  
    633       CALL init_obfbdata( fbdata ) 
    634  
    635       CALL alloc_obfbdata( fbdata, 1, seaicedata%nsurf, 1, 1, 0, .TRUE. ) 
    636  
    637       fbdata%cname(1)      = 'SEAICE' 
    638       fbdata%coblong(1)    = 'Sea ice' 
    639       fbdata%cobunit(1)    = 'Fraction' 
    640       DO je = 1, next 
    641          fbdata%cextname(je) = pext%cdname(je) 
    642          fbdata%cextlong(je) = pext%cdlong(je,1) 
    643          fbdata%cextunit(je) = pext%cdunit(je,1) 
    644       END DO 
    645       fbdata%caddname(1)   = 'Hx' 
    646       fbdata%caddlong(1,1) = 'Model interpolated ICE' 
    647       fbdata%caddunit(1,1) = 'Fraction' 
    648       fbdata%cgrid(1)      = 'T' 
    649       DO ja = 1, nadd 
    650          fbdata%caddname(1+ja) = padd%cdname(ja) 
    651          fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    652          fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    653       END DO 
    654  
    655       WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
    656  
    657       IF(lwp) THEN 
    658          WRITE(numout,*) 
    659          WRITE(numout,*)'obs_wri_seaice :' 
    660          WRITE(numout,*)'~~~~~~~~~~~~~~~~' 
    661          WRITE(numout,*)'Writing SEAICE feedback file : ',TRIM(cfname) 
    662       ENDIF 
    663  
    664       ! Transform obs_prof data structure into obfbdata structure 
    665       fbdata%cdjuldref = '19500101000000' 
    666       DO jo = 1, seaicedata%nsurf 
    667          fbdata%plam(jo)      = seaicedata%rlam(jo) 
    668          fbdata%pphi(jo)      = seaicedata%rphi(jo) 
    669          WRITE(fbdata%cdtyp(jo),'(I4)') seaicedata%ntyp(jo) 
    670          fbdata%ivqc(jo,:)    = 0 
    671          fbdata%ivqcf(:,jo,:) = 0 
    672          IF ( seaicedata%nqc(jo) > 10 ) THEN 
    673             fbdata%ioqc(jo)    = 4 
    674             fbdata%ioqcf(1,jo) = 0 
    675             fbdata%ioqcf(2,jo) = seaicedata%nqc(jo) - 10 
    676          ELSE 
    677             fbdata%ioqc(jo)    = MAX(seaicedata%nqc(jo),1) 
    678             fbdata%ioqcf(:,jo) = 0 
    679          ENDIF 
    680          fbdata%ipqc(jo)      = 0 
    681          fbdata%ipqcf(:,jo)   = 0 
    682          fbdata%itqc(jo)      = 0 
    683          fbdata%itqcf(:,jo)   = 0 
    684          fbdata%cdwmo(jo)     = '' 
    685          fbdata%kindex(jo)    = seaicedata%nsfil(jo) 
    686          IF (ln_grid_global) THEN 
    687             fbdata%iobsi(jo,1) = seaicedata%mi(jo) 
    688             fbdata%iobsj(jo,1) = seaicedata%mj(jo) 
    689          ELSE 
    690             fbdata%iobsi(jo,1) = mig(seaicedata%mi(jo)) 
    691             fbdata%iobsj(jo,1) = mjg(seaicedata%mj(jo)) 
    692          ENDIF 
    693          CALL greg2jul( 0, & 
    694             &           seaicedata%nmin(jo), & 
    695             &           seaicedata%nhou(jo), & 
    696             &           seaicedata%nday(jo), & 
    697             &           seaicedata%nmon(jo), & 
    698             &           seaicedata%nyea(jo), & 
    699             &           fbdata%ptim(jo),   & 
    700             &           krefdate = 19500101 ) 
    701          fbdata%padd(1,jo,1,1) = seaicedata%rmod(jo,1) 
    702          fbdata%pob(1,jo,1)    = seaicedata%robs(jo,1) 
    703          fbdata%pdep(1,jo)     = 0.0 
    704          fbdata%idqc(1,jo)     = 0 
    705          fbdata%idqcf(:,1,jo)  = 0 
    706          IF ( seaicedata%nqc(jo) > 10 ) THEN 
    707             fbdata%ivlqc(1,jo,1) = 4 
    708             fbdata%ivlqcf(1,1,jo,1) = 0 
    709             fbdata%ivlqcf(2,1,jo,1) = seaicedata%nqc(jo) - 10 
    710          ELSE 
    711             fbdata%ivlqc(1,jo,1) = MAX(seaicedata%nqc(jo),1) 
    712             fbdata%ivlqcf(:,1,jo,1) = 0 
    713          ENDIF 
    714          fbdata%iobsk(1,jo,1)  = 0 
    715          DO ja = 1, nadd 
    716             fbdata%padd(1,jo,1+ja,1) = & 
    717                & seaicedata%rext(jo,padd%ipoint(ja)) 
    718          END DO 
    719          DO je = 1, next 
    720             fbdata%pext(1,jo,je) = & 
    721                & seaicedata%rext(jo,pext%ipoint(je)) 
    722          END DO 
    723  
    724       END DO 
    725  
    726       ! Write the obfbdata structure 
    727       CALL write_obfbdata( cfname, fbdata ) 
    728  
    729       ! Output some basic statistics 
    730       CALL obs_wri_stats( fbdata ) 
    731  
    732       CALL dealloc_obfbdata( fbdata ) 
    733  
    734    END SUBROUTINE obs_wri_seaice 
    735  
    736    SUBROUTINE obs_wri_vel( cprefix, profdata, k2dint, padd, pext ) 
    737       !!----------------------------------------------------------------------- 
    738       !! 
    739       !!                     *** ROUTINE obs_wri_vel  *** 
    740       !! 
    741       !! ** Purpose : Write current (profile) observation  
    742       !!              related diagnostics 
    743       !! 
    744       !! ** Method  : NetCDF 
    745       !!  
    746       !! ** Action  : 
    747       !! 
    748       !! History : 
    749       !!      ! 09-01  (K. Mogensen) New feedback format routine 
    750       !!----------------------------------------------------------------------- 
    751  
    752       !! * Modules used 
    753  
    754       !! * Arguments 
    755       CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
    756       TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data 
    757       INTEGER, INTENT(IN) :: k2dint                 ! Horizontal interpolation method 
    758       TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
    759       TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
    760  
    761       !! * Local declarations 
    762       TYPE(obfbdata) :: fbdata 
    763       CHARACTER(LEN=40) :: cfname 
    764       INTEGER :: ilevel 
    765       INTEGER :: jvar 
    766       INTEGER :: jk 
    767       INTEGER :: ik 
    768       INTEGER :: jo 
    769       INTEGER :: ja 
    770       INTEGER :: je 
    771       INTEGER :: nadd 
    772       INTEGER :: next 
    773       REAL(wp) :: zpres 
    774       REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
    775          & zu, & 
    776          & zv 
    777  
    778       IF ( PRESENT( padd ) ) THEN 
    779          nadd = padd%inum 
    780       ELSE 
    781          nadd = 0 
    782       ENDIF 
    783  
    784       IF ( PRESENT( pext ) ) THEN 
    785          next = pext%inum 
    786       ELSE 
    787          next = 0 
    788       ENDIF 
    789  
    790       CALL init_obfbdata( fbdata ) 
    791  
    792       ! Find maximum level 
    793       ilevel = 0 
    794       DO jvar = 1, 2 
    795          ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 
    796       END DO 
    797       CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. ) 
    798  
    799       fbdata%cname(1)      = 'UVEL' 
    800       fbdata%cname(2)      = 'VVEL' 
    801       fbdata%coblong(1)    = 'Zonal velocity' 
    802       fbdata%coblong(2)    = 'Meridional velocity' 
    803       fbdata%cobunit(1)    = 'm/s' 
    804       fbdata%cobunit(2)    = 'm/s' 
    805       DO je = 1, next 
    806          fbdata%cextname(je) = pext%cdname(je) 
    807          fbdata%cextlong(je) = pext%cdlong(je,1) 
    808          fbdata%cextunit(je) = pext%cdunit(je,1) 
    809       END DO 
    810       fbdata%caddname(1)   = 'Hx' 
    811       fbdata%caddlong(1,1) = 'Model interpolated zonal velocity' 
    812       fbdata%caddlong(1,2) = 'Model interpolated meridional velocity' 
    813       fbdata%caddunit(1,1) = 'm/s' 
    814       fbdata%caddunit(1,2) = 'm/s' 
    815       fbdata%caddname(2)   = 'HxG' 
    816       fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)' 
    817       fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)' 
    818       fbdata%caddunit(2,1) = 'm/s' 
    819       fbdata%caddunit(2,2) = 'm/s'  
    820       fbdata%cgrid(1)      = 'U'  
    821       fbdata%cgrid(2)      = 'V' 
    822       DO ja = 1, nadd 
    823          fbdata%caddname(2+ja) = padd%cdname(ja) 
    824          fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1) 
    825          fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1) 
    826       END DO 
    827  
    828       WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
    829  
    830       IF(lwp) THEN 
    831          WRITE(numout,*) 
    832          WRITE(numout,*)'obs_wri_vel :' 
    833          WRITE(numout,*)'~~~~~~~~~~~~~' 
    834          WRITE(numout,*)'Writing velocuty feedback file : ',TRIM(cfname) 
    835       ENDIF 
    836  
    837       ALLOCATE( & 
    838          & zu(profdata%nvprot(1)), & 
    839          & zv(profdata%nvprot(2))  & 
    840          & ) 
    841       CALL obs_rotvel( profdata, k2dint, zu, zv ) 
    842  
    843       ! Transform obs_prof data structure into obfbdata structure 
    844       fbdata%cdjuldref = '19500101000000' 
    845       DO jo = 1, profdata%nprof 
    846          fbdata%plam(jo)      = profdata%rlam(jo) 
    847          fbdata%pphi(jo)      = profdata%rphi(jo) 
    848          WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo) 
    849          fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:) 
    850          fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 
    851          IF ( profdata%nqc(jo) > 10 ) THEN 
    852             fbdata%ioqc(jo)    = 4 
    853             fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 
    854             fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10 
    855          ELSE 
    856             fbdata%ioqc(jo)    = profdata%nqc(jo) 
    857             fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo) 
    858          ENDIF 
    859          fbdata%ipqc(jo)      = profdata%ipqc(jo) 
    860          fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo) 
    861          fbdata%itqc(jo)      = profdata%itqc(jo) 
    862          fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo) 
    863          fbdata%cdwmo(jo)     = profdata%cwmo(jo) 
    864          fbdata%kindex(jo)    = profdata%npfil(jo) 
    865          DO jvar = 1, profdata%nvar 
    866             IF (ln_grid_global) THEN 
    867                fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar) 
    868                fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar) 
    869             ELSE 
    870                fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar)) 
    871                fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar)) 
    872             ENDIF 
    873          END DO 
    874          CALL greg2jul( 0, & 
    875             &           profdata%nmin(jo), & 
    876             &           profdata%nhou(jo), & 
    877             &           profdata%nday(jo), & 
    878             &           profdata%nmon(jo), & 
    879             &           profdata%nyea(jo), & 
    880             &           fbdata%ptim(jo),   & 
    881             &           krefdate = 19500101 ) 
    882          ! Reform the profiles arrays for output 
    883          DO jvar = 1, 2 
    884             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    885                ik = profdata%var(jvar)%nvlidx(jk) 
    886                IF ( jvar == 1 ) THEN 
    887                   fbdata%padd(ik,jo,1,jvar) = zu(jk) 
    888                ELSE 
    889                   fbdata%padd(ik,jo,1,jvar) = zv(jk) 
    890                ENDIF 
    891                fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk) 
    892                fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk) 
    893                fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk) 
    894                fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk) 
    895                fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk) 
    896                IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN 
    897                   fbdata%ivlqc(ik,jo,jvar) = 4 
    898                   fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 
    899                   fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 
    900                ELSE 
    901                   fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) 
    902                   fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk) 
    903                ENDIF 
    904                fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk) 
    905                DO ja = 1, nadd 
    906                   fbdata%padd(ik,jo,2+ja,jvar) = & 
    907                      & profdata%var(jvar)%vext(jk,padd%ipoint(ja)) 
    908                END DO 
    909                DO je = 1, next 
    910                   fbdata%pext(ik,jo,je) = & 
    911                      & profdata%var(jvar)%vext(jk,pext%ipoint(je)) 
    912                END DO 
    913             END DO 
    914          END DO 
    915       END DO 
    916  
    917       ! Write the obfbdata structure 
    918       CALL write_obfbdata( cfname, fbdata ) 
    919        
    920       ! Output some basic statistics 
    921       CALL obs_wri_stats( fbdata ) 
    922  
    923       CALL dealloc_obfbdata( fbdata ) 
    924       
    925       DEALLOCATE( & 
    926          & zu, & 
    927          & zv  & 
    928          & ) 
    929  
    930    END SUBROUTINE obs_wri_vel 
     535   END SUBROUTINE obs_wri_surf 
    931536 
    932537   SUBROUTINE obs_wri_stats( fbdata ) 
     
    952557      INTEGER :: jk 
    953558 
    954 !      INTEGER :: nlev 
    955 !      INTEGER :: nlevmpp 
    956 !      INTEGER :: nobsmpp 
    957559      INTEGER :: numgoodobs 
    958560      INTEGER :: numgoodobsmpp 
Note: See TracChangeset for help on using the changeset viewer.