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 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 – NEMO

Ignore:
Timestamp:
2015-12-21T12:35:23+01:00 (8 years ago)
Author:
timgraham
Message:

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r4990 r6140  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   'key_diaobs' : Switch on the observation diagnostic computation 
    10    !!---------------------------------------------------------------------- 
    119   !!   dia_obs_init : Reading and prepare observations 
    1210   !!   dia_obs      : Compute model equivalent to observations 
    1311   !!   dia_obs_wri  : Write observational diagnostics 
     12   !!   calc_date    : Compute the date of timestep in YYYYMMDD.HHMMSS format 
    1413   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS 
    1514   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS 
    1615   !!---------------------------------------------------------------------- 
    17    !! * Modules used    
     16   !! * Modules used 
    1817   USE wrk_nemo                 ! Memory Allocation 
    1918   USE par_kind                 ! Precision variables 
    … …  
    2120   USE par_oce 
    2221   USE dom_oce                  ! Ocean space and time domain variables 
    23    USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch 
    24    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   
     22   USE obs_read_prof            ! Reading and allocation of profile obs 
     23   USE obs_read_surf            ! Reading and allocation of surface obs 
     24   USE obs_sstbias              ! Bias correction routine for SST  
    2725   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 
    3026   USE obs_prep                 ! Preparation of obs. (grid search etc). 
    3127   USE obs_oper                 ! Observation operators 
    … …  
    3430   USE obs_read_altbias         ! Bias treatment for altimeter 
    3531   USE obs_profiles_def         ! Profile data definitions 
    36    USE obs_profiles             ! Profile data storage 
    3732   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 
    4133   USE obs_types                ! Definitions for observation types 
    4234   USE mpp_map                  ! MPP mapping 
    … …  
    5042      &   dia_obs,      &  ! Compute model equivalent to observations 
    5143      &   dia_obs_wri,  &  ! Write model equivalent to observations 
    52       &   dia_obs_dealloc  ! Deallocate dia_obs data 
    53  
    54    !! * Shared Module variables 
    55    LOGICAL, PUBLIC, PARAMETER :: & 
    56 #if defined key_diaobs 
    57       & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics 
    58 #else 
    59       & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics 
    60 #endif 
     44      &   dia_obs_dealloc, &  ! Deallocate dia_obs data 
     45      &   calc_date           ! Compute the date of a timestep 
    6146 
    6247   !! * Module variables 
    63    LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles 
    64    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 
    68    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 
    71    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 
    75    LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration 
    76    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 
    82    LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height 
    83    LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity 
    84    LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations 
    85    LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land 
    86    LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias   
    87    LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files 
    88    LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations 
    89  
    90    REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS 
    91    REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS 
    92    
    93    INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method 
    94    INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method  
    95  
     48   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator 
     49   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs 
     50    
     51   INTEGER :: nn_1dint       !: Vertical interpolation method 
     52   INTEGER :: nn_2dint       !: Horizontal interpolation method 
    9653   INTEGER, DIMENSION(imaxavtypes) :: & 
    97       & endailyavtypes !: ENACT data types which are daily average 
    98  
    99    INTEGER, PARAMETER :: MaxNumFiles = 1000 
    100    LOGICAL, DIMENSION(MaxNumFiles) :: & 
    101       & ln_profb_ena, & !: Is the feedback files from ENACT data ? 
    102    !                    !: If so use endailyavtypes 
    103       & ln_profb_enatim !: Change tim for 820 enact data set. 
    104     
    105    LOGICAL, DIMENSION(MaxNumFiles) :: & 
    106       & ln_velfb_av   !: Is the velocity feedback files daily average? 
    107    LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    108       & ld_enact     !: Profile data is ENACT so use endailyavtypes 
    109    LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    110       & ld_velav     !: Velocity data is daily averaged 
    111    LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    112       & ld_sstnight  !: SST observation corresponds to night mean 
     54      & nn_profdavtypes      !: Profile data types representing a daily average 
     55   INTEGER :: nproftypes     !: Number of profile obs types 
     56   INTEGER :: nsurftypes     !: Number of surface obs types 
     57   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     58      & nvarsprof, &         !: Number of profile variables 
     59      & nvarssurf            !: Number of surface variables 
     60   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     61      & nextrprof, &         !: Number of profile extra variables 
     62      & nextrsurf            !: Number of surface extra variables 
     63   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: sstbias_type !SST bias type     
     64   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
     65      & surfdata, &          !: Initial surface data 
     66      & surfdataqc           !: Surface data after quality control 
     67   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 
     68      & profdata, &          !: Initial profile data 
     69      & profdataqc           !: Profile data after quality control 
     70 
     71   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 
     72      & cobstypesprof, &     !: Profile obs types 
     73      & cobstypessurf        !: Surface obs types 
    11374 
    11475   !!---------------------------------------------------------------------- 
    … …  
    13596      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    13697      !!        !  07-03  (K. Mogensen) General handling of profiles 
     98      !!        !  14-08  (J.While) Incorporated SST bias correction   
     99      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    137100      !!---------------------------------------------------------------------- 
    138101 
    … …  
    140103 
    141104      !! * Local declarations 
    142       CHARACTER(len=128) :: enactfiles(MaxNumFiles) 
    143       CHARACTER(len=128) :: coriofiles(MaxNumFiles) 
    144       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)       
    149       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) 
    157       CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 
    158       CHARACTER(LEN=128) :: reysstname 
    159       CHARACTER(LEN=12)  :: reysstfmt 
    160       CHARACTER(LEN=128) :: bias_file 
    161       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,                         & 
    169          &            dobsini, dobsend, n1dint, n2dint,               & 
    170          &            nmsshc, mdtcorr, mdtcutoff,                     & 
    171          &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    & 
    172          &            ln_sstnight,                                    & 
     105      INTEGER, PARAMETER :: & 
     106         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type 
     107      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     108         & ifilesprof, &         ! Number of profile files 
     109         & ifilessurf            ! Number of surface files 
     110      INTEGER :: ios             ! Local integer output status for namelist read 
     111      INTEGER :: jtype           ! Counter for obs types 
     112      INTEGER :: jvar            ! Counter for variables 
     113      INTEGER :: jfile           ! Counter for files 
     114 
     115      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
     116         & cn_profbfiles, &      ! T/S profile input filenames 
     117         & cn_sstfbfiles, &      ! Sea surface temperature input filenames 
     118         & cn_slafbfiles, &      ! Sea level anomaly input filenames 
     119         & cn_sicfbfiles, &      ! Seaice concentration input filenames 
     120         & cn_velfbfiles, &      ! Velocity profile input filenames 
     121         & cn_sstbias_files      ! SST bias input filenames 
     122      CHARACTER(LEN=128) :: & 
     123         & cn_altbiasfile        ! Altimeter bias input filename 
     124      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
     125         & clproffiles, &        ! Profile filenames 
     126         & clsurffiles           ! Surface filenames 
     127 
     128      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
     129      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
     130      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
     131      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
     132      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     133      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     134      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
     135      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
     136      LOGICAL :: ln_sstbias     !: Logical switch for bias corection of SST  
     137      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
     138      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     139      LOGICAL :: llvar1          ! Logical for profile variable 1 
     140      LOGICAL :: llvar2          ! Logical for profile variable 1 
     141      LOGICAL :: llnightav       ! Logical for calculating night-time averages 
     142      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
     143 
     144      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
     145      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
     146      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     147         & zglam1, &             ! Model longitudes for profile variable 1 
     148         & zglam2                ! Model longitudes for profile variable 2 
     149      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     150         & zgphi1, &             ! Model latitudes for profile variable 1 
     151         & zgphi2                ! Model latitudes for profile variable 2 
     152      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     153         & zmask1, &             ! Model land/sea mask associated with variable 1 
     154         & zmask2                ! Model land/sea mask associated with variable 2 
     155 
     156      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
     157         &            ln_sst, ln_sic, ln_vel3d,                       & 
     158         &            ln_altbias, ln_nea, ln_grid_global,             & 
    173159         &            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 
    203       INTEGER :: ios                 ! Local integer output status for namelist read 
    204       LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 
     160         &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     161         &            cn_profbfiles, cn_slafbfiles,                   & 
     162         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
     163         &            cn_velfbfiles, cn_altbiasfile,                  & 
     164         &            cn_gridsearchfile, rn_gridsearchres,            & 
     165         &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     166         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
     167         &            nn_profdavtypes, ln_sstbias, cn_sstbias_files 
     168 
     169      INTEGER :: jnumsstbias 
     170      CALL wrk_alloc( jpi, jpj, zglam1 ) 
     171      CALL wrk_alloc( jpi, jpj, zglam2 ) 
     172      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
     173      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
     174      CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 
     175      CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 
    205176 
    206177      !----------------------------------------------------------------------- 
    207178      ! Read namelist parameters 
    208179      !----------------------------------------------------------------------- 
    209  
    210       enactfiles(:) = '' 
    211       coriofiles(:) = '' 
    212       profbfiles(:) = '' 
    213       slafilesact(:) = '' 
    214       slafilespas(:) = '' 
    215       slafbfiles(:) = '' 
    216       sstfiles(:)   = '' 
    217       sstfbfiles(:) = '' 
    218       seaicefiles(:) = '' 
    219       velcurfiles(:) = '' 
    220       veladcpfiles(:) = '' 
    221       velavcurfiles(:) = '' 
    222       velhrcurfiles(:) = '' 
    223       velavadcpfiles(:) = '' 
    224       velhradcpfiles(:) = '' 
    225       velfbfiles(:) = '' 
    226       velcurfiles(:) = '' 
    227       veladcpfiles(:) = '' 
    228       endailyavtypes(:) = -1 
    229       endailyavtypes(1) = 820 
    230       ln_profb_ena(:) = .FALSE. 
    231       ln_profb_enatim(:) = .TRUE. 
    232       ln_velfb_av(:) = .FALSE. 
    233       ln_ignmis = .FALSE. 
    234180       
    235       CALL ini_date( dobsini ) 
    236       CALL fin_date( dobsend ) 
    237   
    238       ! Read Namelist namobs : control observation diagnostics 
    239       REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation 
     181      !Initalise all values in namelist arrays 
     182      ALLOCATE(sstbias_type(jpmaxnfiles)) 
     183      ! Some namelist arrays need initialising 
     184      cn_profbfiles(:) = '' 
     185      cn_slafbfiles(:) = '' 
     186      cn_sstfbfiles(:) = '' 
     187      cn_sicfbfiles(:) = '' 
     188      cn_velfbfiles(:) = '' 
     189      cn_sstbias_files(:) = '' 
     190      nn_profdavtypes(:) = -1 
     191 
     192      CALL ini_date( rn_dobsini ) 
     193      CALL fin_date( rn_dobsend ) 
     194 
     195      ! Read namelist namobs : control observation diagnostics 
     196      REWIND( numnam_ref )   ! Namelist namobs in reference namelist 
    240197      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 
    241198901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 
    242199 
    243       REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation 
     200      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist 
    244201      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 
    245202902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 
    246203      IF(lwm) WRITE ( numond, namobs ) 
    247204 
    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) 
     205      IF ( .NOT. ln_diaobs ) THEN 
     206         IF(lwp) WRITE(numout,cform_war) 
     207         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs' 
     208         RETURN 
     209      ENDIF 
     210       
     211      !----------------------------------------------------------------------- 
     212      ! Set up list of observation types to be used 
     213      ! and the files associated with each type 
     214      !----------------------------------------------------------------------- 
     215 
     216      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     217      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
     218 
     219      IF (ln_sstbias) THEN  
     220         lmask(:) = .FALSE.  
     221         WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
     222         jnumsstbias = COUNT(lmask)  
     223         lmask(:) = .FALSE.  
    282224      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 
     225 
     226      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     227         IF(lwp) WRITE(numout,cform_war) 
     228         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     229            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
     230            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     231         nwarn = nwarn + 1 
     232         ln_diaobs = .FALSE. 
     233         RETURN 
     234      ENDIF 
     235 
     236      IF ( nproftypes > 0 ) THEN 
     237 
     238         ALLOCATE( cobstypesprof(nproftypes) ) 
     239         ALLOCATE( ifilesprof(nproftypes) ) 
     240         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     241 
     242         jtype = 0 
     243         IF (ln_t3d .OR. ln_s3d) THEN 
     244            jtype = jtype + 1 
     245            clproffiles(jtype,:) = cn_profbfiles(:) 
     246            cobstypesprof(jtype) = 'prof  ' 
     247            ifilesprof(jtype) = 0 
     248            DO jfile = 1, jpmaxnfiles 
     249               IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
     250                  ifilesprof(jtype) = ifilesprof(jtype) + 1 
     251            END DO 
     252         ENDIF 
     253         IF (ln_vel3d) THEN 
     254            jtype = jtype + 1 
     255            clproffiles(jtype,:) = cn_velfbfiles(:) 
     256            cobstypesprof(jtype) = 'vel   ' 
     257            ifilesprof(jtype) = 0 
     258            DO jfile = 1, jpmaxnfiles 
     259               IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
     260                  ifilesprof(jtype) = ifilesprof(jtype) + 1 
     261            END DO 
     262         ENDIF 
     263 
     264      ENDIF 
     265 
     266      IF ( nsurftypes > 0 ) THEN 
     267 
     268         ALLOCATE( cobstypessurf(nsurftypes) ) 
     269         ALLOCATE( ifilessurf(nsurftypes) ) 
     270         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     271 
     272         jtype = 0 
     273         IF (ln_sla) THEN 
     274            jtype = jtype + 1 
     275            clsurffiles(jtype,:) = cn_slafbfiles(:) 
     276            cobstypessurf(jtype) = 'sla   ' 
     277            ifilessurf(jtype) = 0 
     278            DO jfile = 1, jpmaxnfiles 
     279               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     280                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     281            END DO 
     282         ENDIF 
     283         IF (ln_sst) THEN 
     284            jtype = jtype + 1 
     285            clsurffiles(jtype,:) = cn_sstfbfiles(:) 
     286            cobstypessurf(jtype) = 'sst   ' 
     287            ifilessurf(jtype) = 0 
     288            DO jfile = 1, jpmaxnfiles 
     289               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     290                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     291            END DO 
     292         ENDIF 
     293#if defined key_lim2 || defined key_lim3 
     294         IF (ln_sic) THEN 
     295            jtype = jtype + 1 
     296            clsurffiles(jtype,:) = cn_sicfbfiles(:) 
     297            cobstypessurf(jtype) = 'sic   ' 
     298            ifilessurf(jtype) = 0 
     299            DO jfile = 1, jpmaxnfiles 
     300               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     301                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     302            END DO 
     303         ENDIF 
     304#endif 
     305 
     306      ENDIF 
     307 
     308      !Write namelist settings to stdout 
    322309      IF(lwp) THEN 
    323310         WRITE(numout,*) 
    … …  
    325312         WRITE(numout,*) '~~~~~~~~~~~~' 
    326313         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters'  
    327          WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d 
    328          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 
    332          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 
    335          WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh 
    336          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 
    340          WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight 
    341          WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss 
    342          WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice 
    343          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 
    350          WRITE(numout,*) & 
    351    '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup 
     314         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d 
     315         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d 
     316         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
     317         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
     318         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
     319         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
     320         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
     321         WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias  
     322         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
    352323         IF (ln_grid_search_lookup) & 
    353             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)) 
     324            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     325         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini 
     326         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
     327         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
     328         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
     329         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     330         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
     331         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
     332         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
     333         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     334         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
     335         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
     336         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
     337         WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     338 
     339         IF ( nproftypes > 0 ) THEN 
     340            DO jtype = 1, nproftypes 
     341               DO jfile = 1, ifilesprof(jtype) 
     342                  WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
     343                     TRIM(clproffiles(jtype,jfile)) 
     344               END DO 
    358345            END DO 
    359346         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)) 
     347 
     348         WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
     349         IF ( nsurftypes > 0 ) THEN 
     350            DO jtype = 1, nsurftypes 
     351               DO jfile = 1, ifilessurf(jtype) 
     352                  WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
     353                     TRIM(clsurffiles(jtype,jfile)) 
     354               END DO 
    364355            END DO 
    365356         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 
    448          WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend 
    449          WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint 
    450          WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint 
    451          WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea 
    452          WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc 
    453          WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr 
    454          WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff 
    455          WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias 
    456          WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis 
    457          WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes 
    458  
    459       ENDIF 
    460        
     357         WRITE(numout,*) '~~~~~~~~~~~~' 
     358 
     359      ENDIF 
     360 
     361      !----------------------------------------------------------------------- 
     362      ! Obs operator parameter checking and initialisations 
     363      !----------------------------------------------------------------------- 
     364 
    461365      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 
    462366         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) 
    … …  
    464368      ENDIF 
    465369 
    466       CALL obs_typ_init 
    467        
    468       CALL mppmap_init 
    469        
    470       ! Parameter control 
    471 #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 
    476          IF(lwp) WRITE(numout,cform_war) 
    477          IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 
    478             &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 
    479          nwarn = nwarn + 1 
    480       ENDIF 
    481 #endif 
    482  
    483       CALL obs_grid_setup( ) 
    484       IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 
     370      IF ( ln_grid_global ) THEN 
     371         CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' ) 
     372      ENDIF 
     373 
     374      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 
    485375         CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 
    486376            &                    ' is not available') 
    487377      ENDIF 
    488       IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 
     378 
     379      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
    489380         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    490381            &                    ' is not available') 
    491382      ENDIF 
    492383 
     384      CALL obs_typ_init 
     385      IF(ln_grid_global) THEN 
     386         CALL mppmap_init 
     387      ENDIF 
     388 
     389      CALL obs_grid_setup( ) 
     390 
    493391      !----------------------------------------------------------------------- 
    494392      ! Depending on switches read the various observation types 
    495393      !----------------------------------------------------------------------- 
    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 
    524              
    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  
     394 
     395      IF ( nproftypes > 0 ) THEN 
     396 
     397         ALLOCATE(profdata(nproftypes)) 
     398         ALLOCATE(profdataqc(nproftypes)) 
     399         ALLOCATE(nvarsprof(nproftypes)) 
     400         ALLOCATE(nextrprof(nproftypes)) 
     401 
     402         DO jtype = 1, nproftypes 
     403 
     404            nvarsprof(jtype) = 2 
     405            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     406               nextrprof(jtype) = 1 
     407               llvar1 = ln_t3d 
     408               llvar2 = ln_s3d 
     409               zglam1 = glamt 
     410               zgphi1 = gphit 
     411               zmask1 = tmask 
     412               zglam2 = glamt 
     413               zgphi2 = gphit 
     414               zmask2 = tmask 
     415            ENDIF 
     416            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     417               nextrprof(jtype) = 2 
     418               llvar1 = ln_vel3d 
     419               llvar2 = ln_vel3d 
     420               zglam1 = glamu 
     421               zgphi1 = gphiu 
     422               zmask1 = umask 
     423               zglam2 = glamv 
     424               zgphi2 = gphiv 
     425               zmask2 = vmask 
     426            ENDIF 
     427 
     428            !Read in profile or profile obs types 
     429            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       & 
     430               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
     431               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
     432               &               rn_dobsini, rn_dobsend, llvar1, llvar2, & 
     433               &               ln_ignmis, ln_s_at_t, .FALSE., & 
     434               &               kdailyavtypes = nn_profdavtypes ) 
     435 
     436            DO jvar = 1, nvarsprof(jtype) 
     437               CALL obs_prof_staend( profdata(jtype), jvar ) 
    539438            END DO 
    540439 
    541             CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    542                &              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 
     440            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
     441               &               llvar1, llvar2, & 
     442               &               jpi, jpj, jpk, & 
     443               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
     444               &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     445 
     446         END DO 
     447 
     448         DEALLOCATE( ifilesprof, clproffiles ) 
     449 
     450      ENDIF 
     451 
     452      IF ( nsurftypes > 0 ) THEN 
     453 
     454         ALLOCATE(surfdata(nsurftypes)) 
     455         ALLOCATE(surfdataqc(nsurftypes)) 
     456         ALLOCATE(nvarssurf(nsurftypes)) 
     457         ALLOCATE(nextrsurf(nsurftypes)) 
     458 
     459         DO jtype = 1, nsurftypes 
     460 
     461            nvarssurf(jtype) = 1 
     462            nextrsurf(jtype) = 0 
     463            llnightav = .FALSE. 
     464            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
     465            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 
     466 
     467            !Read in surface obs types 
     468            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
     469               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
     470               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
     471               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    620472          
    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 
    629473          
    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       
     474            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     475 
     476            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     477               CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
     478               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     479            ENDIF 
     480            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
     481               !Read in bias field and correct SST. 
     482               IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
     483                                                     "  but no bias"// & 
     484                                                     " files to read in")    
     485                  CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
     486                                        jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 
     487            ENDIF 
     488         END DO 
     489 
     490         DEALLOCATE( ifilessurf, clsurffiles ) 
     491 
     492      ENDIF 
     493 
     494      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
     495      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
     496      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
     497      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
     498      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 
     499      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 
     500 
    967501   END SUBROUTINE dia_obs_init 
    968502 
    … …  
    974508      !! 
    975509      !! ** Method  : Call the observation operators on each time step to 
    976       !!              compute the model equivalent of the following date: 
    977       !!               - T profiles 
    978       !!               - S profiles 
    979       !!               - Sea surface height (referenced to a mean) 
    980       !!               - Sea surface temperature 
    981       !!               - Sea surface salinity 
    982       !!               - Velocity component (U,V) profiles 
    983       !! 
    984       !! ** Action  :  
     510      !!              compute the model equivalent of the following data: 
     511      !!               - Profile data, currently T/S or U/V 
     512      !!               - Surface data, currently SST, SLA or sea-ice concentration. 
     513      !! 
     514      !! ** Action  : 
    985515      !! 
    986516      !! History : 
    … …  
    991521      !!        !  07-04  (G. Smith) Generalized surface operators 
    992522      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles 
     523      !!        !  14-08  (J. While) observation operator for profiles in  
     524      !!                             generalised vertical coordinates 
     525      !!        !  15-08  (M. Martin) Combined surface/profile routines. 
    993526      !!---------------------------------------------------------------------- 
    994527      !! * Modules used 
    995528      USE dom_oce, ONLY : &             ! Ocean space and time domain variables 
    996          & rdt,           &                        
    997          & gdept_1d,       &              
    998          & tmask, umask, vmask                             
     529         & gdept_n,       &       
     530         & gdept_1d       
    999531      USE phycst, ONLY : &              ! Physical constants 
    1000532         & rday                          
    1001533      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables 
    1002534         & tsn,  &              
    1003          & un, vn,  & 
    1004          & sshn 
     535         & un, vn, & 
     536         & sshn   
     537      USE phycst, ONLY : &         ! Physical constants 
     538         & rday 
    1005539#if defined  key_lim3 
    1006       USE ice, ONLY : &                     ! LIM Ice model variables 
     540      USE ice, ONLY : &            ! LIM3 Ice model variables 
    1007541         & frld 
    1008542#endif 
    1009543#if defined key_lim2 
    1010       USE ice_2, ONLY : &                     ! LIM Ice model variables 
     544      USE ice_2, ONLY : &          ! LIM2 Ice model variables 
    1011545         & frld 
    1012546#endif 
    … …  
    1014548 
    1015549      !! * Arguments 
    1016       INTEGER, INTENT(IN) :: kstp                         ! Current timestep 
     550      INTEGER, INTENT(IN) :: kstp  ! Current timestep 
    1017551      !! * Local declarations 
    1018       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 
    1024       INTEGER :: jvar                   ! Variable number     
     552      INTEGER :: idaystp           ! Number of timesteps per day 
     553      INTEGER :: jtype             ! Data loop variable 
     554      INTEGER :: jvar              ! Variable number 
     555      INTEGER :: ji, jj            ! Loop counters 
     556      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     557         & zprofvar1, &            ! Model values for 1st variable in a prof ob 
     558         & zprofvar2               ! Model values for 2nd variable in a prof ob 
     559      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     560         & zprofmask1, &           ! Mask associated with zprofvar1 
     561         & zprofmask2              ! Mask associated with zprofvar2 
     562      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     563         & zsurfvar                ! Model values equivalent to surface ob. 
     564      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     565         & zglam1,    &            ! Model longitudes for prof variable 1 
     566         & zglam2,    &            ! Model longitudes for prof variable 2 
     567         & zgphi1,    &            ! Model latitudes for prof variable 1 
     568         & zgphi2                  ! Model latitudes for prof variable 2 
    1025569#if ! defined key_lim2 && ! defined key_lim3 
    1026       REAL(wp), POINTER, DIMENSION(:,:) :: frld    
     570      REAL(wp), POINTER, DIMENSION(:,:) :: frld 
    1027571#endif 
    1028       CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    1029   
     572      LOGICAL :: llnightav        ! Logical for calculating night-time average 
     573 
     574      !Allocate local work arrays 
     575      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 
     576      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 
     577      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 
     578      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
     579      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     580      CALL wrk_alloc( jpi, jpj, zglam1 ) 
     581      CALL wrk_alloc( jpi, jpj, zglam2 ) 
     582      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
     583      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    1030584#if ! defined key_lim2 && ! defined key_lim3 
    1031585      CALL wrk_alloc(jpi,jpj,frld)  
    … …  
    1047601#endif 
    1048602      !----------------------------------------------------------------------- 
    1049       ! Depending on switches call various observation operators 
    1050       !----------------------------------------------------------------------- 
    1051  
    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 ) 
     603      ! Call the profile and surface observation operators 
     604      !----------------------------------------------------------------------- 
     605 
     606      IF ( nproftypes > 0 ) THEN 
     607 
     608         DO jtype = 1, nproftypes 
     609 
     610            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
     611            CASE('prof') 
     612               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 
     613               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 
     614               zprofmask1(:,:,:) = tmask(:,:,:) 
     615               zprofmask2(:,:,:) = tmask(:,:,:) 
     616               zglam1(:,:) = glamt(:,:) 
     617               zglam2(:,:) = glamt(:,:) 
     618               zgphi1(:,:) = gphit(:,:) 
     619               zgphi2(:,:) = gphit(:,:) 
     620            CASE('vel') 
     621               zprofvar1(:,:,:) = un(:,:,:) 
     622               zprofvar2(:,:,:) = vn(:,:,:) 
     623               zprofmask1(:,:,:) = umask(:,:,:) 
     624               zprofmask2(:,:,:) = vmask(:,:,:) 
     625               zglam1(:,:) = glamu(:,:) 
     626               zglam2(:,:) = glamv(:,:) 
     627               zgphi1(:,:) = gphiu(:,:) 
     628               zgphi2(:,:) = gphiv(:,:) 
     629            END SELECT 
     630 
     631            IF( ln_zco .OR. ln_zps ) THEN  
     632               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     633                  &               nit000, idaystp,                         & 
     634                  &               zprofvar1, zprofvar2,                    & 
     635                  &               gdept_1d, zprofmask1, zprofmask2,        & 
     636                  &               zglam1, zglam2, zgphi1, zgphi2,          & 
     637                  &               nn_1dint, nn_2dint,                      & 
     638                  &               kdailyavtypes = nn_profdavtypes ) 
     639            ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 
     640               !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 
     641               CALL obs_pro_sco_opt( profdataqc(jtype),                    &  
     642                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
     643                  &              zprofvar1, zprofvar2,                   &  
     644                  &              gdept_n(:,:,:), gdepw_n(:,:,:),           & 
     645                  &              tmask, nn_1dint, nn_2dint,              &  
     646                  &              kdailyavtypes = nn_profdavtypes )  
    1061647            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              ) 
     648               CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 
     649                         'yet working for velocity data (turn off velocity observations') 
    1066650            ENDIF 
     651 
    1067652         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) ) 
     653 
     654      ENDIF 
     655 
     656      IF ( nsurftypes > 0 ) THEN 
     657 
     658         DO jtype = 1, nsurftypes 
     659 
     660            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
     661            CASE('sst') 
     662               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
     663               llnightav = ln_sstnight 
     664            CASE('sla') 
     665               zsurfvar(:,:) = sshn(:,:) 
     666               llnightav = .FALSE. 
     667#if defined key_lim2 || defined key_lim3 
     668            CASE('sic') 
     669               IF ( kstp == 0 ) THEN 
     670                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
     671                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
     672                        &           'time-step but some obs are valid then.' ) 
     673                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
     674                        &           ' sea-ice obs will be missed' 
     675                  ENDIF 
     676                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
     677                     &                        surfdataqc(jtype)%nsstp(1) 
     678                  CYCLE 
     679               ELSE 
     680                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     681               ENDIF 
     682 
     683               llnightav = .FALSE. 
     684#endif 
     685            END SELECT 
     686 
     687            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     688               &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
     689               &               nn_2dint, llnightav ) 
     690 
    1086691         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       
     692 
     693      ENDIF 
     694 
     695      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 
     696      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 
     697      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 
     698      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
     699      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     700      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
     701      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
     702      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
     703      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
     704#if ! defined key_lim2 && ! defined key_lim3 
     705      CALL wrk_dealloc(jpi,jpj,frld) 
    1102706#endif 
    1103707 
    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  
    1114 #if ! defined key_lim2 && ! defined key_lim3 
    1115       CALL wrk_dealloc(jpi,jpj,frld)  
    1116 #endif 
    1117  
    1118708   END SUBROUTINE dia_obs 
    1119    
    1120    SUBROUTINE dia_obs_wri  
     709 
     710   SUBROUTINE dia_obs_wri 
    1121711      !!---------------------------------------------------------------------- 
    1122712      !!                    ***  ROUTINE dia_obs_wri  *** 
    … …  
    1126716      !! ** Method  : Call observation diagnostic output routines 
    1127717      !! 
    1128       !! ** Action  :  
     718      !! ** Action  : 
    1129719      !! 
    1130720      !! History : 
    … …  
    1134724      !!        !  07-03  (K. Mogensen) General handling of profiles 
    1135725      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles 
    1136       !!---------------------------------------------------------------------- 
     726      !!        !  15-08  (M. Martin) Combined writing for prof and surf types 
     727      !!---------------------------------------------------------------------- 
     728      !! * Modules used 
     729      USE obs_rot_vel          ! Rotation of velocities 
     730 
    1137731      IMPLICIT NONE 
    1138732 
    1139733      !! * Local declarations 
    1140  
    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 
     734      INTEGER :: jtype                    ! Data set loop variable 
     735      INTEGER :: jo, jvar, jk 
     736      REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     737         & zu, & 
     738         & zv 
     739 
    1150740      !----------------------------------------------------------------------- 
    1151741      ! Depending on switches call various observation output routines 
    1152742      !----------------------------------------------------------------------- 
    1153743 
    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 ) 
     744      IF ( nproftypes > 0 ) THEN 
     745 
     746         DO jtype = 1, nproftypes 
     747 
     748            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 
     749 
     750               ! For velocity data, rotate the model velocities to N/S, E/W 
     751               ! using the compressed data structure. 
     752               ALLOCATE( & 
     753                  & zu(profdataqc(jtype)%nvprot(1)), & 
     754                  & zv(profdataqc(jtype)%nvprot(2))  & 
     755                  & ) 
     756 
     757               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 
     758 
     759               DO jo = 1, profdataqc(jtype)%nprof 
     760                  DO jvar = 1, 2 
     761                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 
     762 
     763                        IF ( jvar == 1 ) THEN 
     764                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 
     765                        ELSE 
     766                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 
     767                        ENDIF 
     768 
     769                     END DO 
     770                  END DO 
     771               END DO 
     772 
     773               DEALLOCATE( zu ) 
     774               DEALLOCATE( zv ) 
     775 
     776            END IF 
     777 
     778            CALL obs_prof_decompress( profdataqc(jtype), & 
     779               &                      profdata(jtype), .TRUE., numout ) 
     780 
     781            CALL obs_wri_prof( profdata(jtype) ) 
    1163782 
    1164783         END DO 
    1165784 
    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 ) 
     785      ENDIF 
     786 
     787      IF ( nsurftypes > 0 ) THEN 
     788 
     789         DO jtype = 1, nsurftypes 
     790 
     791            CALL obs_surf_decompress( surfdataqc(jtype), & 
     792               &                      surfdata(jtype), .TRUE., numout ) 
     793 
     794            CALL obs_wri_surf( surfdata(jtype) ) 
    1216795 
    1217796         END DO 
    1218797 
    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           
    1390798      ENDIF 
    1391799 
    … …  
    1405813      !! 
    1406814      !!---------------------------------------------------------------------- 
    1407       !! obs_grid deallocation 
     815      ! obs_grid deallocation 
    1408816      CALL obs_grid_deallocate 
    1409817 
    1410       !! 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 
     818      ! diaobs deallocation 
     819      IF ( nproftypes > 0 ) & 
     820         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 
     821 
     822      IF ( nsurftypes > 0 ) & 
     823         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     824 
    1433825   END SUBROUTINE dia_obs_dealloc 
    1434826 
    1435    SUBROUTINE ini_date( ddobsini ) 
    1436       !!---------------------------------------------------------------------- 
    1437       !!                    ***  ROUTINE ini_date  *** 
     827   SUBROUTINE calc_date( kstp, ddobs ) 
     828      !!---------------------------------------------------------------------- 
     829      !!                    ***  ROUTINE calc_date  *** 
    1438830      !!           
    1439       !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format 
    1440       !! 
    1441       !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format 
    1442       !! 
    1443       !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format 
     831      !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format 
     832      !! 
     833      !! ** Method  : Get date in double precision YYYYMMDD.HHMMSS format 
     834      !! 
     835      !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format 
     836      !! 
     837      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format 
    1444838      !! 
    1445839      !! History : 
    … …  
    1449843      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date 
    1450844      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2 
     845      !!        !  2014-09  (D. Lea) New generic routine now deals with arbitrary initial time of day 
    1451846      !!---------------------------------------------------------------------- 
    1452847      USE phycst, ONLY : &            ! Physical constants 
    1453848         & rday 
    1454 !      USE daymod, ONLY : &            ! Time variables 
    1455 !         & nmonth_len            
    1456849      USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    1457850         & rdt 
    … …  
    1460853 
    1461854      !! * Arguments 
    1462       REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS 
     855      REAL(KIND=dp), INTENT(OUT) :: ddobs                        ! Date in YYYYMMDD.HHMMSS 
     856      INTEGER :: kstp 
    1463857 
    1464858      !! * Local declarations 
    … …  
    1468862      INTEGER :: ihou 
    1469863      INTEGER :: imin 
    1470       INTEGER :: imday         ! Number of days in month. 
    1471       REAL(KIND=wp) :: zdayfrc ! Fraction of day 
     864      INTEGER :: imday       ! Number of days in month. 
     865      REAL(wp) :: zdayfrc    ! Fraction of day 
    1472866 
    1473867      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year 
    … …  
    1475869      !!---------------------------------------------------------------------- 
    1476870      !! Initial date initialization (year, month, day, hour, minute) 
    1477       !! (This assumes that the initial date is for 00z)) 
    1478871      !!---------------------------------------------------------------------- 
    1479872      iyea =   ndate0 / 10000 
    1480873      imon = ( ndate0 - iyea * 10000 ) / 100 
    1481874      iday =   ndate0 - iyea * 10000 - imon * 100 
    1482       ihou = 0 
    1483       imin = 0 
     875      ihou =   nn_time0 / 100 
     876      imin = ( nn_time0 - ihou * 100 )  
    1484877 
    1485878      !!---------------------------------------------------------------------- 
    1486879      !! Compute number of days + number of hours + min since initial time 
    1487880      !!---------------------------------------------------------------------- 
    1488       iday = iday + ( nit000 -1 ) * rdt / rday 
    1489       zdayfrc = ( nit000 -1 ) * rdt / rday 
     881      zdayfrc = kstp * rdt / rday 
    1490882      zdayfrc = zdayfrc - aint(zdayfrc) 
    1491       ihou = int( zdayfrc * 24 ) 
    1492       imin = int( (zdayfrc * 24 - ihou) * 60 ) 
    1493  
    1494       !!----------------------------------------------------------------------- 
    1495       !! Convert number of days (iday) into a real date 
    1496       !!---------------------------------------------------------------------- 
     883      imin = imin + int( zdayfrc * 24 * 60 )  
     884      DO WHILE (imin >= 60)  
     885        imin=imin-60 
     886        ihou=ihou+1 
     887      END DO 
     888      DO WHILE (ihou >= 24) 
     889        ihou=ihou-24 
     890        iday=iday+1 
     891      END DO  
     892      iday = iday + kstp * rdt / rday  
     893 
     894      !----------------------------------------------------------------------- 
     895      ! Convert number of days (iday) into a real date 
     896      !---------------------------------------------------------------------- 
    1497897 
    1498898      CALL calc_month_len( iyea, imonth_len ) 
    1499        
     899 
    1500900      DO WHILE ( iday > imonth_len(imon) ) 
    1501901         iday = iday - imonth_len(imon) 
    … …  
    1508908      END DO 
    1509909 
    1510       !!---------------------------------------------------------------------- 
    1511       !! Convert it into YYYYMMDD.HHMMSS format. 
    1512       !!---------------------------------------------------------------------- 
    1513       ddobsini = iyea * 10000_dp + imon * 100_dp + & 
    1514          &       iday + ihou * 0.01_dp + imin * 0.0001_dp 
    1515  
    1516  
    1517    END SUBROUTINE ini_date 
    1518  
    1519    SUBROUTINE fin_date( ddobsfin ) 
    1520       !!---------------------------------------------------------------------- 
    1521       !!                    ***  ROUTINE fin_date  *** 
     910      !---------------------------------------------------------------------- 
     911      ! Convert it into YYYYMMDD.HHMMSS format. 
     912      !---------------------------------------------------------------------- 
     913      ddobs = iyea * 10000_dp + imon * 100_dp + & 
     914         &    iday + ihou * 0.01_dp + imin * 0.0001_dp 
     915 
     916   END SUBROUTINE calc_date 
     917 
     918   SUBROUTINE ini_date( ddobsini ) 
     919      !!---------------------------------------------------------------------- 
     920      !!                    ***  ROUTINE ini_date  *** 
    1522921      !!           
    1523       !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format 
    1524       !! 
    1525       !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format 
    1526       !! 
    1527       !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format 
     922      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 
     923      !! 
     924      !! ** Method  :  
     925      !! 
     926      !! ** Action  :  
    1528927      !! 
    1529928      !! History : 
    … …  
    1532931      !!        !  06-10  (A. Weaver) Cleaning 
    1533932      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2 
    1534       !!---------------------------------------------------------------------- 
    1535       USE phycst, ONLY : &            ! Physical constants 
    1536          & rday 
    1537 !      USE daymod, ONLY : &            ! Time variables 
    1538 !         & nmonth_len                 
    1539       USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    1540          & rdt 
     933      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date 
     934      !!---------------------------------------------------------------------- 
    1541935 
    1542936      IMPLICIT NONE 
    1543937 
    1544938      !! * Arguments 
    1545       REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS 
    1546  
    1547       !! * Local declarations 
    1548       INTEGER :: iyea        ! date - (year, month, day, hour, minute) 
    1549       INTEGER :: imon 
    1550       INTEGER :: iday 
    1551       INTEGER :: ihou 
    1552       INTEGER :: imin 
    1553       INTEGER :: imday         ! Number of days in month. 
    1554       REAL(KIND=wp) :: zdayfrc       ! Fraction of day 
    1555           
    1556       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year 
    1557              
    1558       !----------------------------------------------------------------------- 
    1559       ! Initial date initialization (year, month, day, hour, minute) 
    1560       ! (This assumes that the initial date is for 00z) 
    1561       !----------------------------------------------------------------------- 
    1562       iyea =   ndate0 / 10000 
    1563       imon = ( ndate0 - iyea * 10000 ) / 100 
    1564       iday =   ndate0 - iyea * 10000 - imon * 100 
    1565       ihou = 0 
    1566       imin = 0 
    1567        
    1568       !----------------------------------------------------------------------- 
    1569       ! Compute number of days + number of hours + min since initial time 
    1570       !----------------------------------------------------------------------- 
    1571       iday    = iday +  nitend  * rdt / rday 
    1572       zdayfrc =  nitend  * rdt / rday 
    1573       zdayfrc = zdayfrc - AINT( zdayfrc ) 
    1574       ihou    = INT( zdayfrc * 24 ) 
    1575       imin    = INT( ( zdayfrc * 24 - ihou ) * 60 ) 
    1576  
    1577       !----------------------------------------------------------------------- 
    1578       ! Convert number of days (iday) into a real date 
    1579       !---------------------------------------------------------------------- 
    1580  
    1581       CALL calc_month_len( iyea, imonth_len ) 
    1582        
    1583       DO WHILE ( iday > imonth_len(imon) ) 
    1584          iday = iday - imonth_len(imon) 
    1585          imon = imon + 1  
    1586          IF ( imon > 12 ) THEN 
    1587             imon = 1 
    1588             iyea = iyea + 1 
    1589             CALL calc_month_len( iyea, imonth_len )  ! update month lengths 
    1590          ENDIF 
    1591       END DO 
    1592  
    1593       !----------------------------------------------------------------------- 
    1594       ! Convert it into YYYYMMDD.HHMMSS format 
    1595       !----------------------------------------------------------------------- 
    1596       ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday & 
    1597          &     + ihou * 0.01_dp  + imin * 0.0001_dp 
    1598  
    1599     END SUBROUTINE fin_date 
    1600      
     939      REAL(KIND=dp), INTENT(OUT) :: ddobsini                   ! Initial date in YYYYMMDD.HHMMSS 
     940 
     941      CALL calc_date( nit000 - 1, ddobsini ) 
     942 
     943   END SUBROUTINE ini_date 
     944 
     945   SUBROUTINE fin_date( ddobsfin ) 
     946      !!---------------------------------------------------------------------- 
     947      !!                    ***  ROUTINE fin_date  *** 
     948      !!           
     949      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 
     950      !! 
     951      !! ** Method  :  
     952      !! 
     953      !! ** Action  :  
     954      !! 
     955      !! History : 
     956      !!        !  06-03  (K. Mogensen)  Original code 
     957      !!        !  06-05  (K. Mogensen)  Reformatted 
     958      !!        !  06-10  (A. Weaver) Cleaning 
     959      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2 
     960      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date 
     961      !!---------------------------------------------------------------------- 
     962 
     963      IMPLICIT NONE 
     964 
     965      !! * Arguments 
     966      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 
     967 
     968      CALL calc_date( nitend, ddobsfin ) 
     969 
     970   END SUBROUTINE fin_date 
     971    
    1601972END MODULE diaobs 
Note: See TracChangeset for help on using the changeset viewer.