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

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

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

Updated OBS simplification branch to the head of the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r4245 r5659  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   obs_pro_opt :    Compute the model counterpart of temperature and 
    10    !!                    salinity observations from profiles 
    11    !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly 
    12    !!                    observations 
    13    !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature 
    14    !!                    observations 
    15    !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity 
    16    !!                    observations 
    17    !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration 
    18    !!                    observations 
    19    !! 
    20    !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional 
    21    !!                    components of velocity from observations. 
     9   !!   obs_prof_opt :    Compute the model counterpart of profile data 
     10   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    2211   !!---------------------------------------------------------------------- 
    2312 
    … …  
    4635   PRIVATE 
    4736 
    48    PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations 
    49       &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations 
    50       &   obs_sst_opt, &  ! Compute the model counterpart of SST observations 
    51       &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations 
    52       &   obs_seaice_opt, & 
    53       &   obs_vel_opt     ! Compute the model counterpart of velocity profile data 
     37   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile observations 
     38      &   obs_surf_opt     ! Compute the model counterpart of surface observations 
    5439 
    5540   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 
    … …  
    6348CONTAINS 
    6449 
    65    SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    66       &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, & 
    67       &                    kdailyavtypes ) 
     50   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
     51      &                     pvar1, pvar2, pgdept, pmask1, k1dint, k2dint, & 
     52      &                     kdailyavtypes ) 
    6853      !!----------------------------------------------------------------------- 
    6954      !! 
    … …  
    7863      !! 
    7964      !!    First, a vertical profile of horizontally interpolated model 
    80       !!    now temperatures is computed at the obs (lon, lat) point. 
     65      !!    now values is computed at the obs (lon, lat) point. 
    8166      !!    Several horizontal interpolation schemes are available: 
    8267      !!        - distance-weighted (great circle) (k2dint = 0) 
    … …  
    8671      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    8772      !! 
    88       !!    Next, the vertical temperature profile is interpolated to the 
     73      !!    Next, the vertical profile is interpolated to the 
    8974      !!    data depth points. Two vertical interpolation schemes are 
    9075      !!    available: 
    … …  
    9681      !!    routine. 
    9782      !! 
    98       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is 
     83      !!    If the logical is switched on, the model equivalent is 
    9984      !!    a daily mean model temperature field. So, we first compute 
    10085      !!    the mean, then interpolate only at the end of the day. 
    10186      !! 
    102       !!    Note: the in situ temperature observations must be converted 
     87      !!    Note: in situ temperature observations must be converted 
    10388      !!    to potential temperature (the model variable) prior to 
    10489      !!    assimilation.  
    105       !!?????????????????????????????????????????????????????????????? 
    106       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 
    107       !!?????????????????????????????????????????????????????????????? 
    10890      !! 
    10991      !! ** Action  : 
    … …  
    11597      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity 
    11698      !!      ! 07-03 (K. Mogensen) General handling of profiles 
     99      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
    117100      !!----------------------------------------------------------------------- 
    118101   
    … …  
    134117      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
    135118      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    136          & ptn,    &    ! Model temperature field 
    137          & psn,    &    ! Model salinity field 
    138          & ptmask       ! Land-sea mask 
     119         & pvar1,    &    ! Model field 1 
     120         & pvar2,    &    ! Model field 2 
     121         & pmask1,   &    ! Land-sea mask 1 
     122         & pmask2         ! Land-sea mask 2 
    139123      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    140124         & pgdept       ! Model array of depth levels 
    … …  
    158142      REAL(KIND=wp) :: zdaystp 
    159143      REAL(KIND=wp), DIMENSION(kpk) :: & 
     144         & zobsmask1, & 
     145         & zobsmask2, & 
    160146         & zobsmask, & 
    161147         & zobsk,    & 
    162148         & zobs2k 
    163149      REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
    164          & zweig 
     150         & zweig1, & 
     151         & zweig2 
    165152      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    166          & zmask, & 
    167          & zintt, & 
    168          & zints, & 
    169          & zinmt, & 
    170          & zinms 
     153         & zmask1, & 
     154         & zmask2, & 
     155         & zint1, & 
     156         & zint2, & 
     157         & zinm1, & 
     158         & zinm2 
    171159      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    172          & zglam, & 
    173          & zgphi 
     160         & zglam1, & 
     161         & zglam2, & 
     162         & zgphi1, & 
     163         & zgphi2 
    174164      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    175          & igrdi, & 
    176          & igrdj 
     165         & igrdi1, & 
     166         & igrdi2, & 
     167         & igrdj1, & 
     168         & igrdj2 
    177169 
    178170      !------------------------------------------------------------------------ 
    … …  
    209201         DO jj = 1, jpj 
    210202            DO ji = 1, jpi 
    211                ! Increment the temperature field for computing daily mean 
     203               ! Increment field 1 for computing daily mean 
    212204               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    213                   &                        + ptn(ji,jj,jk) 
    214                ! Increment the salinity field for computing daily mean 
     205                  &                        + pvar1(ji,jj,jk) 
     206               ! Increment field 2 for computing daily mean 
    215207               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    216                   &                        + psn(ji,jj,jk) 
     208                  &                        + pvar2(ji,jj,jk) 
    217209            END DO 
    218210         END DO 
    … …  
    236228      ! Get the data for interpolation 
    237229      ALLOCATE( & 
    238          & igrdi(2,2,ipro),      & 
    239          & igrdj(2,2,ipro),      & 
    240          & zglam(2,2,ipro),      & 
    241          & zgphi(2,2,ipro),      & 
    242          & zmask(2,2,kpk,ipro),  & 
    243          & zintt(2,2,kpk,ipro),  & 
    244          & zints(2,2,kpk,ipro)   & 
     230         & igrdi1(2,2,ipro),      & 
     231         & igrdi2(2,2,ipro),      & 
     232         & igrdj1(2,2,ipro),      & 
     233         & igrdj2(2,2,ipro),      & 
     234         & zglam1(2,2,ipro),      & 
     235         & zglam2(2,2,ipro),      & 
     236         & zgphi1(2,2,ipro),      & 
     237         & zgphi2(2,2,ipro),      & 
     238         & zmask1(2,2,kpk,ipro),  & 
     239         & zmask2(2,2,kpk,ipro),  & 
     240         & zint1(2,2,kpk,ipro),  & 
     241         & zint2(2,2,kpk,ipro)   & 
    245242         & ) 
    246243 
    247244      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    248245         iobs = jobs - prodatqc%nprofup 
    249          igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    250          igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    251          igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    252          igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 
    253          igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 
    254          igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    255          igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 
    256          igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 
     246         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 
     247         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 
     248         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 
     249         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 
     250         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 
     251         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 
     252         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
     253         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
     254         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
     255         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
     256         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
     257         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
     258         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
     259         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
     260         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
     261         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
    257262      END DO 
    258263 
    259       CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 
    260       CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 
    261       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 
    262       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt ) 
    263       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints ) 
     264      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, glamt1, zglam1 ) 
     265      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, gphit1, zgphi1 ) 
     266      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 
     267      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
     268       
     269      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, glamt2, zglam2 ) 
     270      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, gphit2, zgphi2 ) 
     271      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask1, zmask2 ) 
     272      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    264273 
    265274      ! At the end of the day also get interpolated means 
    … …  
    267276 
    268277         ALLOCATE( & 
    269             & zinmt(2,2,kpk,ipro),  & 
    270             & zinms(2,2,kpk,ipro)   & 
     278            & zinm1(2,2,kpk,ipro),  & 
     279            & zinm2(2,2,kpk,ipro)   & 
    271280            & ) 
    272281 
    273          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    274             &                  prodatqc%vdmean(:,:,:,1), zinmt ) 
    275          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    276             &                  prodatqc%vdmean(:,:,:,2), zinms ) 
     282         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, & 
     283            &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
     284         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, & 
     285            &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
    277286 
    278287      ENDIF 
    … …  
    304313         ! Horizontal weights and vertical mask 
    305314 
    306          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
    307             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 
     315         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    308316 
    309317            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    310                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    311                &                   zmask(:,:,:,iobs), zweig, zobsmask ) 
    312  
     318               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
     319               &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     320 
     321         ENDIF 
     322 
     323         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     324 
     325            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     326               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
     327               &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     328  
    313329         ENDIF 
    314330 
    … …  
    321337               IF ( idayend == 0 )  THEN 
    322338                   
    323                   ! Daily averaged moored buoy (MRB) data 
     339                  ! Daily averaged data 
    324340                   
    325341                  CALL obs_int_h2d( kpk, kpk,      & 
    326                      &              zweig, zinmt(:,:,:,iobs), zobsk ) 
     342                     &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    327343                   
    328344                   
    … …  
    330346                
    331347                  CALL ctl_stop( ' A nonzero' //     & 
    332                      &           ' number of profile T BUOY data should' // & 
     348                     &           ' number of average profile data should' // & 
    333349                     &           ' only occur at the end of a given day' ) 
    334350 
    … …  
    340356 
    341357               CALL obs_int_h2d( kpk, kpk,      & 
    342                   &              zweig, zintt(:,:,:,iobs), zobsk ) 
     358                  &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    343359 
    344360            ENDIF 
    … …  
    377393               IF ( idayend == 0 )  THEN 
    378394 
    379                   ! Daily averaged moored buoy (MRB) data 
     395                  ! Daily averaged data 
    380396                   
    381397                  CALL obs_int_h2d( kpk, kpk,      & 
    382                      &              zweig, zinms(:,:,:,iobs), zobsk ) 
     398                     &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    383399                   
    384400               ELSE 
    385401 
    386402                  CALL ctl_stop( ' A nonzero' //     & 
    387                      &           ' number of profile S BUOY data should' // & 
     403                     &           ' number of average profile data should' // & 
    388404                     &           ' only occur at the end of a given day' ) 
    389405 
    … …  
    395411 
    396412               CALL obs_int_h2d( kpk, kpk,      & 
    397                   &              zweig, zints(:,:,:,iobs), zobsk ) 
     413                  &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    398414 
    399415            ENDIF 
    … …  
    429445      ! Deallocate the data for interpolation 
    430446      DEALLOCATE( & 
    431          & igrdi, & 
    432          & igrdj, & 
    433          & zglam, & 
    434          & zgphi, & 
    435          & zmask, & 
    436          & zintt, & 
    437          & zints  & 
     447         & igrdi1, & 
     448         & igrdi2, & 
     449         & igrdj1, & 
     450         & igrdj2, & 
     451         & zglam1, & 
     452         & zglam2, & 
     453         & zgphi1, & 
     454         & zgphi2, & 
     455         & zmask1, & 
     456         & zmask2, & 
     457         & zint1,  & 
     458         & zint2   & 
    438459         & ) 
    439460      ! At the end of the day also get interpolated means 
    440461      IF ( idayend == 0 ) THEN 
    441462         DEALLOCATE( & 
    442             & zinmt,  & 
    443             & zinms   & 
     463            & zinm1,  & 
     464            & zinm2   & 
    444465            & ) 
    445466      ENDIF 
    … …  
    447468      prodatqc%nprofup = prodatqc%nprofup + ipro  
    448469       
    449    END SUBROUTINE obs_pro_opt 
    450  
    451    SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 
    452       &                    psshn, psshmask, k2dint ) 
     470   END SUBROUTINE obs_prof_opt 
     471 
     472   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, kit000, & 
     473      &                    psurf, psurfmask, k2dint, ld_nightav ) 
    453474      !!----------------------------------------------------------------------- 
    454475      !! 
    455       !!                     ***  ROUTINE obs_sla_opt  *** 
    456       !! 
    457       !! ** Purpose : Compute the model counterpart of sea level anomaly 
     476      !!                     ***  ROUTINE obs_surf_opt  *** 
     477      !! 
     478      !! ** Purpose : Compute the model counterpart of surface 
    458479      !!              data by interpolating from the model grid to the  
    459480      !!              observation point. 
    … …  
    462483      !!              the model values at the corners of the surrounding grid box. 
    463484      !! 
    464       !!    The now model SSH is first computed at the obs (lon, lat) point. 
     485      !!    The new model value is first computed at the obs (lon, lat) point. 
    465486      !! 
    466487      !!    Several horizontal interpolation schemes are available: 
    … …  
    471492      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    472493      !!   
    473       !!    The sea level anomaly at the observation points is then computed  
    474       !!    by removing a mean dynamic topography (defined at the obs. point). 
    475494      !! 
    476495      !! ** Action  : 
    … …  
    478497      !! History : 
    479498      !!      ! 07-03 (A. Weaver) 
     499      !!      ! 15-02 (M. Martin) Combined routine for surface types 
    480500      !!----------------------------------------------------------------------- 
    481501   
    … …  
    486506 
    487507      !! * Arguments 
    488       TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening 
     508      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc     ! Subset of surface data not failing screening 
    489509      INTEGER, INTENT(IN) :: kt      ! Time step 
    490510      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
    … …  
    494514      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    495515      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    496          & psshn,  &    ! Model SSH field 
    497          & psshmask     ! Land-sea mask 
     516         & psurf,  &    ! Model surface field 
     517         & psurfmask    ! Land-sea mask 
    498518          
    499519      !! * Local declarations 
    … …  
    502522      INTEGER :: jobs 
    503523      INTEGER :: inrc 
    504       INTEGER :: isla 
     524      INTEGER :: isurf 
    505525      INTEGER :: iobs 
    506526      REAL(KIND=wp) :: zlam 
    … …  
    511531      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    512532         & zmask, & 
    513          & zsshl, & 
     533         & zsurf, & 
    514534         & zglam, & 
    515535         & zgphi 
    … …  
    523543      ! ... Record and data counters 
    524544      inrc = kt - kit000 + 2 
    525       isla = sladatqc%nsstp(inrc) 
     545      isurf = surfdataqc%nsstp(inrc) 
     546       
     547      IF ( ld_nightav ) THEN 
     548 
     549      ! Initialize array for night mean 
     550         IF ( kt .EQ. 0 ) THEN 
     551            ALLOCATE ( icount_night(kpi,kpj) ) 
     552            ALLOCATE ( imask_night(kpi,kpj) ) 
     553            ALLOCATE ( zintmp(kpi,kpj) ) 
     554            ALLOCATE ( zouttmp(kpi,kpj) ) 
     555            ALLOCATE ( zmeanday(kpi,kpj) ) 
     556            nday_qsr = -1   ! initialisation flag for nbc_dcy 
     557         ENDIF 
     558 
     559         ! Initialize daily mean for first timestep 
     560         idayend = MOD( kt - kit000 + 1, kdaystp ) 
     561 
     562         ! Added kt == 0 test to catch restart case  
     563         IF ( idayend == 1 .OR. kt == 0) THEN 
     564            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 
     565            DO jj = 1, jpj 
     566               DO ji = 1, jpi 
     567                  surfdataqc%vdmean(ji,jj) = 0.0 
     568                  zmeanday(ji,jj) = 0.0 
     569                  icount_night(ji,jj) = 0 
     570               END DO 
     571            END DO 
     572         ENDIF 
     573 
     574         zintmp(:,:) = 0.0 
     575         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
     576         imask_night(:,:) = INT( zouttmp(:,:) ) 
     577 
     578         DO jj = 1, jpj 
     579            DO ji = 1, jpi 
     580               ! Increment the temperature field for computing night mean and counter 
     581               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  & 
     582                      &                        + psurf(ji,jj)*imask_night(ji,jj) 
     583               zmeanday(ji,jj)        = zmeanday(ji,jj) + psurf(ji,jj) 
     584               icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj) 
     585            END DO 
     586         END DO 
     587    
     588         ! Compute the daily mean at the end of day 
     589 
     590         zdaystp = 1.0 / REAL( kdaystp ) 
     591 
     592         IF ( idayend == 0 ) THEN  
     593            DO jj = 1, jpj 
     594               DO ji = 1, jpi 
     595                  ! Test if "no night" point 
     596                  IF ( icount_night(ji,jj) .NE. 0 ) THEN 
     597                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     598                       &                        / icount_night(ji,jj)  
     599                  ELSE 
     600                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
     601                  ENDIF 
     602               END DO 
     603            END DO 
     604         ENDIF 
     605 
     606      ENDIF 
    526607 
    527608      ! Get the data for interpolation 
    528609 
    529610      ALLOCATE( & 
    530          & igrdi(2,2,isla), & 
    531          & igrdj(2,2,isla), & 
    532          & zglam(2,2,isla), & 
    533          & zgphi(2,2,isla), & 
    534          & zmask(2,2,isla), & 
    535          & zsshl(2,2,isla)  & 
     611         & igrdi(2,2,isurf), & 
     612         & igrdj(2,2,isurf), & 
     613         & zglam(2,2,isurf), & 
     614         & zgphi(2,2,isurf), & 
     615         & zmask(2,2,isurf), & 
     616         & zsurf(2,2,isurf)  & 
    536617         & ) 
    537618       
    538       DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 
    539          iobs = jobs - sladatqc%nsurfup 
    540          igrdi(1,1,iobs) = sladatqc%mi(jobs)-1 
    541          igrdj(1,1,iobs) = sladatqc%mj(jobs)-1 
    542          igrdi(1,2,iobs) = sladatqc%mi(jobs)-1 
    543          igrdj(1,2,iobs) = sladatqc%mj(jobs) 
    544          igrdi(2,1,iobs) = sladatqc%mi(jobs) 
    545          igrdj(2,1,iobs) = sladatqc%mj(jobs)-1 
    546          igrdi(2,2,iobs) = sladatqc%mi(jobs) 
    547          igrdj(2,2,iobs) = sladatqc%mj(jobs) 
     619      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     620         iobs = jobs - surfdataqc%nsurfup 
     621         igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
     622         igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
     623         igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
     624         igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
     625         igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
     626         igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
     627         igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
     628         igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
    548629      END DO 
    549630 
    550       CALL obs_int_comm_2d( 2, 2, isla, & 
     631      CALL obs_int_comm_2d( 2, 2, isurf, & 
    551632         &                  igrdi, igrdj, glamt, zglam ) 
    552       CALL obs_int_comm_2d( 2, 2, isla, & 
     633      CALL obs_int_comm_2d( 2, 2, isurf, & 
    553634         &                  igrdi, igrdj, gphit, zgphi ) 
    554       CALL obs_int_comm_2d( 2, 2, isla, & 
    555          &                  igrdi, igrdj, psshmask, zmask ) 
    556       CALL obs_int_comm_2d( 2, 2, isla, & 
    557          &                  igrdi, igrdj, psshn, zsshl ) 
     635      CALL obs_int_comm_2d( 2, 2, isurf, & 
     636         &                  igrdi, igrdj, psurfmask, zmask ) 
     637      CALL obs_int_comm_2d( 2, 2, isurf, & 
     638         &                  igrdi, igrdj, psurf, zsurf ) 
     639 
     640      ! At the end of the day get interpolated means 
     641      IF ( idayend == 0 .AND. ld_nightav ) THEN 
     642 
     643         ALLOCATE( & 
     644            & zsurfm(2,2,isurf)  & 
     645            & ) 
     646 
     647         CALL obs_int_comm_2d( 2, 2, isurf, igrdi, igrdj, & 
     648            &               surfdataqc%vdmean(:,:), zsurfm ) 
     649 
     650      ENDIF 
    558651 
    559652      ! Loop over observations 
    560  
    561       DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 
    562  
    563          iobs = jobs - sladatqc%nsurfup 
    564  
    565          IF ( kt /= sladatqc%mstp(jobs) ) THEN 
     653      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     654 
     655         iobs = jobs - surfdataqc%nsurfup 
     656 
     657         IF ( kt /= surfdataqc%mstp(jobs) ) THEN 
    566658             
    567659            IF(lwp) THEN 
    … …  
    574666               WRITE(numout,*) ' Record  = ', jobs,                & 
    575667                  &            ' kt      = ', kt,                  & 
    576                   &            ' mstp    = ', sladatqc%mstp(jobs), & 
    577                   &            ' ntyp    = ', sladatqc%ntyp(jobs) 
     668                  &            ' mstp    = ', surfdataqc%mstp(jobs), & 
     669                  &            ' ntyp    = ', surfdataqc%ntyp(jobs) 
    578670            ENDIF 
    579671            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 
    … …  
    581673         ENDIF 
    582674          
    583          zlam = sladatqc%rlam(jobs) 
    584          zphi = sladatqc%rphi(jobs) 
    585  
    586          ! Get weights to interpolate the model SSH to the observation point 
     675         zlam = surfdataqc%rlam(jobs) 
     676         zphi = surfdataqc%rphi(jobs) 
     677 
     678         ! Get weights to interpolate the model value to the observation point 
    587679         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    588680            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    … …  
    590682          
    591683 
    592          ! Interpolate the model SSH to the observation point 
    593          CALL obs_int_h2d( 1, 1,      & 
    594             &              zweig, zsshl(:,:,iobs),  zext ) 
     684         ! Interpolate the model field to the observation point 
     685         IF ( ld_nightav ) THEN 
     686 
     687           IF ( idayend == 0 )  THEN 
     688               ! Daily averaged data 
     689               CALL obs_int_h2d( 1, 1,      &  
     690                     &              zweig, zsurfm(:,:,iobs), zext ) 
     691            ELSE  
     692               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     & 
     693                     &           ' number of night data should' // & 
     694                     &           ' only occur at the end of a given day' ) 
     695            ENDIF 
     696 
     697         ELSE 
     698 
     699            CALL obs_int_h2d( 1, 1,      & 
     700               &              zweig, zsurf(:,:,iobs),  zext ) 
     701 
     702         ENDIF 
    595703          
    596          sladatqc%rext(jobs,1) = zext(1) 
    597          ! ... Remove the MDT at the observation point 
    598          sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 
     704         IF ( surfdataqc%nextr == 2 ) THEN 
     705            ! ... Remove the MDT from the SSH at the observation point to get the SLA 
     706            surfdataqc%rext(jobs,1) = zext(1) 
     707            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
     708         ELSE 
     709            surfdataqc%rmod(jobs,1) = zext(1) 
     710         ENDIF 
    599711 
    600712      END DO 
    … …  
    607719         & zgphi, & 
    608720         & zmask, & 
    609          & zsshl  & 
    610          & ) 
    611  
    612       sladatqc%nsurfup = sladatqc%nsurfup + isla 
    613  
    614    END SUBROUTINE obs_sla_opt 
    615  
    616    SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 
    617       &                    psstn, psstmask, k2dint, ld_nightav ) 
    618       !!----------------------------------------------------------------------- 
    619       !! 
    620       !!                     ***  ROUTINE obs_sst_opt  *** 
    621       !! 
    622       !! ** Purpose : Compute the model counterpart of surface temperature 
    623       !!              data by interpolating from the model grid to the  
    624       !!              observation point. 
    625       !! 
    626       !! ** Method  : Linearly interpolate to each observation point using  
    627       !!              the model values at the corners of the surrounding grid box. 
    628       !! 
    629       !!    The now model SST is first computed at the obs (lon, lat) point. 
    630       !! 
    631       !!    Several horizontal interpolation schemes are available: 
    632       !!        - distance-weighted (great circle) (k2dint = 0) 
    633       !!        - distance-weighted (small angle)  (k2dint = 1) 
    634       !!        - bilinear (geographical grid)     (k2dint = 2) 
    635       !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    636       !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    637       !! 
    638       !! 
    639       !! ** Action  : 
    640       !! 
    641       !! History : 
    642       !!        !  07-07  (S. Ricci ) : Original 
    643       !!       
    644       !!----------------------------------------------------------------------- 
    645  
    646       !! * Modules used 
    647       USE obs_surf_def  ! Definition of storage space for surface observations 
    648       USE sbcdcy 
    649  
    650       IMPLICIT NONE 
    651  
    652       !! * Arguments 
    653       TYPE(obs_surf), INTENT(INOUT) :: & 
    654          & sstdatqc     ! Subset of surface data not failing screening 
    655       INTEGER, INTENT(IN) :: kt        ! Time step 
    656       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    657       INTEGER, INTENT(IN) :: kpj 
    658       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    659                                        !   (kit000-1 = restart time) 
    660       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    661       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day   
    662       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    663          & psstn,  &    ! Model SST field 
    664          & psstmask     ! Land-sea mask 
    665  
    666       !! * Local declarations 
    667       INTEGER :: ji 
    668       INTEGER :: jj 
    669       INTEGER :: jobs 
    670       INTEGER :: inrc 
    671       INTEGER :: isst 
    672       INTEGER :: iobs 
    673       INTEGER :: idayend 
    674       REAL(KIND=wp) :: zlam 
    675       REAL(KIND=wp) :: zphi 
    676       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    677       REAL(KIND=wp) :: zdaystp 
    678       INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    679          & icount_sstnight,      & 
    680          & imask_night 
    681       REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    682          & zintmp, & 
    683          & zouttmp, &  
    684          & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    685       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    686          & zweig 
    687       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    688          & zmask, & 
    689          & zsstl, & 
    690          & zsstm, & 
    691          & zglam, & 
    692          & zgphi 
    693       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    694          & igrdi, & 
    695          & igrdj 
    696       LOGICAL, INTENT(IN) :: ld_nightav 
    697  
    698       !----------------------------------------------------------------------- 
    699       ! Local initialization  
    700       !----------------------------------------------------------------------- 
    701       ! ... Record and data counters 
    702       inrc = kt - kit000 + 2 
    703       isst = sstdatqc%nsstp(inrc) 
    704  
    705       IF ( ld_nightav ) THEN 
    706  
    707       ! Initialize array for night mean 
    708  
    709       IF ( kt .EQ. 0 ) THEN 
    710          ALLOCATE ( icount_sstnight(kpi,kpj) ) 
    711          ALLOCATE ( imask_night(kpi,kpj) ) 
    712          ALLOCATE ( zintmp(kpi,kpj) ) 
    713          ALLOCATE ( zouttmp(kpi,kpj) ) 
    714          ALLOCATE ( zmeanday(kpi,kpj) ) 
    715          nday_qsr = -1   ! initialisation flag for nbc_dcy 
    716       ENDIF 
    717  
    718       ! Initialize daily mean for first timestep 
    719       idayend = MOD( kt - kit000 + 1, kdaystp ) 
    720  
    721       ! Added kt == 0 test to catch restart case  
    722       IF ( idayend == 1 .OR. kt == 0) THEN 
    723          IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 
    724          DO jj = 1, jpj 
    725             DO ji = 1, jpi 
    726                sstdatqc%vdmean(ji,jj) = 0.0 
    727                zmeanday(ji,jj) = 0.0 
    728                icount_sstnight(ji,jj) = 0 
    729             END DO 
    730          END DO 
    731       ENDIF 
    732  
    733       zintmp(:,:) = 0.0 
    734       zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
    735       imask_night(:,:) = INT( zouttmp(:,:) ) 
    736  
    737       DO jj = 1, jpj 
    738          DO ji = 1, jpi 
    739             ! Increment the temperature field for computing night mean and counter 
    740             sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  & 
    741                    &                        + psstn(ji,jj)*imask_night(ji,jj) 
    742             zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj) 
    743             icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 
    744          END DO 
    745       END DO 
    746     
    747       ! Compute the daily mean at the end of day 
    748  
    749       zdaystp = 1.0 / REAL( kdaystp ) 
    750  
    751       IF ( idayend == 0 ) THEN  
    752          DO jj = 1, jpj 
    753             DO ji = 1, jpi 
    754                ! Test if "no night" point 
    755                IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 
    756                   sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 
    757                     &                        / icount_sstnight(ji,jj)  
    758                ELSE 
    759                   sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
    760                ENDIF 
    761             END DO 
    762          END DO 
    763       ENDIF 
    764  
    765       ENDIF 
    766  
    767       ! Get the data for interpolation 
    768        
    769       ALLOCATE( & 
    770          & igrdi(2,2,isst), & 
    771          & igrdj(2,2,isst), & 
    772          & zglam(2,2,isst), & 
    773          & zgphi(2,2,isst), & 
    774          & zmask(2,2,isst), & 
    775          & zsstl(2,2,isst)  & 
    776          & ) 
    777        
    778       DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 
    779          iobs = jobs - sstdatqc%nsurfup 
    780          igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 
    781          igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 
    782          igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 
    783          igrdj(1,2,iobs) = sstdatqc%mj(jobs) 
    784          igrdi(2,1,iobs) = sstdatqc%mi(jobs) 
    785          igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 
    786          igrdi(2,2,iobs) = sstdatqc%mi(jobs) 
    787          igrdj(2,2,iobs) = sstdatqc%mj(jobs) 
    788       END DO 
    789        
    790       CALL obs_int_comm_2d( 2, 2, isst, & 
    791          &                  igrdi, igrdj, glamt, zglam ) 
    792       CALL obs_int_comm_2d( 2, 2, isst, & 
    793          &                  igrdi, igrdj, gphit, zgphi ) 
    794       CALL obs_int_comm_2d( 2, 2, isst, & 
    795          &                  igrdi, igrdj, psstmask, zmask ) 
    796       CALL obs_int_comm_2d( 2, 2, isst, & 
    797          &                  igrdi, igrdj, psstn, zsstl ) 
    798  
    799       ! At the end of the day get interpolated means 
    800       IF ( idayend == 0 .AND. ld_nightav ) THEN 
    801  
    802          ALLOCATE( & 
    803             & zsstm(2,2,isst)  & 
    804             & ) 
    805  
    806          CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 
    807             &               sstdatqc%vdmean(:,:), zsstm ) 
    808  
    809       ENDIF 
    810  
    811       ! Loop over observations 
    812  
    813       DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 
    814           
    815          iobs = jobs - sstdatqc%nsurfup 
    816           
    817          IF ( kt /= sstdatqc%mstp(jobs) ) THEN 
    818              
    819             IF(lwp) THEN 
    820                WRITE(numout,*) 
    821                WRITE(numout,*) ' E R R O R : Observation',              & 
    822                   &            ' time step is not consistent with the', & 
    823                   &            ' model time step' 
    824                WRITE(numout,*) ' =========' 
    825                WRITE(numout,*) 
    826                WRITE(numout,*) ' Record  = ', jobs,                & 
    827                   &            ' kt      = ', kt,                  & 
    828                   &            ' mstp    = ', sstdatqc%mstp(jobs), & 
    829                   &            ' ntyp    = ', sstdatqc%ntyp(jobs) 
    830             ENDIF 
    831             CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 
    832              
    833          ENDIF 
    834           
    835          zlam = sstdatqc%rlam(jobs) 
    836          zphi = sstdatqc%rphi(jobs) 
    837           
    838          ! Get weights to interpolate the model SST to the observation point 
    839          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    840             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    841             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    842              
    843          ! Interpolate the model SST to the observation point  
    844  
    845          IF ( ld_nightav ) THEN 
    846  
    847            IF ( idayend == 0 )  THEN 
    848                ! Daily averaged/diurnal cycle of SST  data 
    849                CALL obs_int_h2d( 1, 1,      &  
    850                      &              zweig, zsstm(:,:,iobs), zext ) 
    851             ELSE  
    852                CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     & 
    853                      &           ' number of night SST data should' // & 
    854                      &           ' only occur at the end of a given day' ) 
    855             ENDIF 
    856  
    857          ELSE 
    858  
    859             CALL obs_int_h2d( 1, 1,      & 
    860             &              zweig, zsstl(:,:,iobs),  zext ) 
    861  
    862          ENDIF 
    863          sstdatqc%rmod(jobs,1) = zext(1) 
    864           
    865       END DO 
    866        
    867       ! Deallocate the data for interpolation 
    868       DEALLOCATE( & 
    869          & igrdi, & 
    870          & igrdj, & 
    871          & zglam, & 
    872          & zgphi, & 
    873          & zmask, & 
    874          & zsstl  & 
     721         & zsurf  & 
    875722         & ) 
    876723 
    … …  
    878725      IF ( idayend == 0 .AND. ld_nightav ) THEN 
    879726         DEALLOCATE( & 
    880             & zsstm  & 
     727            & zsurfm  & 
    881728            & ) 
    882729      ENDIF 
    883        
    884       sstdatqc%nsurfup = sstdatqc%nsurfup + isst 
    885  
    886    END SUBROUTINE obs_sst_opt 
    887  
    888    SUBROUTINE obs_sss_opt 
    889       !!----------------------------------------------------------------------- 
    890       !! 
    891       !!                     ***  ROUTINE obs_sss_opt  *** 
    892       !! 
    893       !! ** Purpose : Compute the model counterpart of sea surface salinity 
    894       !!              data by interpolating from the model grid to the  
    895       !!              observation point. 
    896       !! 
    897       !! ** Method  :  
    898       !! 
    899       !! ** Action  : 
    900       !! 
    901       !! History : 
    902       !!      ! ??-??  
    903       !!----------------------------------------------------------------------- 
    904  
    905       IMPLICIT NONE 
    906  
    907    END SUBROUTINE obs_sss_opt 
    908  
    909    SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, & 
    910       &                    pseaicen, pseaicemask, k2dint ) 
    911  
    912       !!----------------------------------------------------------------------- 
    913       !! 
    914       !!                     ***  ROUTINE obs_seaice_opt  *** 
    915       !! 
    916       !! ** Purpose : Compute the model counterpart of surface temperature 
    917       !!              data by interpolating from the model grid to the  
    918       !!              observation point. 
    919       !! 
    920       !! ** Method  : Linearly interpolate to each observation point using  
    921       !!              the model values at the corners of the surrounding grid box. 
    922       !! 
    923       !!    The now model sea ice is first computed at the obs (lon, lat) point. 
    924       !! 
    925       !!    Several horizontal interpolation schemes are available: 
    926       !!        - distance-weighted (great circle) (k2dint = 0) 
    927       !!        - distance-weighted (small angle)  (k2dint = 1) 
    928       !!        - bilinear (geographical grid)     (k2dint = 2) 
    929       !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    930       !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    931       !! 
    932       !! 
    933       !! ** Action  : 
    934       !! 
    935       !! History : 
    936       !!        !  07-07  (S. Ricci ) : Original 
    937       !!       
    938       !!----------------------------------------------------------------------- 
    939  
    940       !! * Modules used 
    941       USE obs_surf_def  ! Definition of storage space for surface observations 
    942  
    943       IMPLICIT NONE 
    944  
    945       !! * Arguments 
    946       TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening 
    947       INTEGER, INTENT(IN) :: kt       ! Time step 
    948       INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
    949       INTEGER, INTENT(IN) :: kpj 
    950       INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
    951                                       !   (kit000-1 = restart time) 
    952       INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    953       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    954          & pseaicen,  &    ! Model sea ice field 
    955          & pseaicemask     ! Land-sea mask 
    956           
    957       !! * Local declarations 
    958       INTEGER :: ji 
    959       INTEGER :: jj 
    960       INTEGER :: jobs 
    961       INTEGER :: inrc 
    962       INTEGER :: iseaice 
    963       INTEGER :: iobs 
    964         
    965       REAL(KIND=wp) :: zlam 
    966       REAL(KIND=wp) :: zphi 
    967       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    968       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    969          & zweig 
    970       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    971          & zmask, & 
    972          & zseaicel, & 
    973          & zglam, & 
    974          & zgphi 
    975       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    976          & igrdi, & 
    977          & igrdj 
    978  
    979       !------------------------------------------------------------------------ 
    980       ! Local initialization  
    981       !------------------------------------------------------------------------ 
    982       ! ... Record and data counters 
    983       inrc = kt - kit000 + 2 
    984       iseaice = seaicedatqc%nsstp(inrc) 
    985  
    986       ! Get the data for interpolation 
    987        
    988       ALLOCATE( & 
    989          & igrdi(2,2,iseaice), & 
    990          & igrdj(2,2,iseaice), & 
    991          & zglam(2,2,iseaice), & 
    992          & zgphi(2,2,iseaice), & 
    993          & zmask(2,2,iseaice), & 
    994          & zseaicel(2,2,iseaice)  & 
    995          & ) 
    996        
    997       DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 
    998          iobs = jobs - seaicedatqc%nsurfup 
    999          igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1 
    1000          igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1 
    1001          igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1 
    1002          igrdj(1,2,iobs) = seaicedatqc%mj(jobs) 
    1003          igrdi(2,1,iobs) = seaicedatqc%mi(jobs) 
    1004          igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1 
    1005          igrdi(2,2,iobs) = seaicedatqc%mi(jobs) 
    1006          igrdj(2,2,iobs) = seaicedatqc%mj(jobs) 
    1007       END DO 
    1008        
    1009       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1010          &                  igrdi, igrdj, glamt, zglam ) 
    1011       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1012          &                  igrdi, igrdj, gphit, zgphi ) 
    1013       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1014          &                  igrdi, igrdj, pseaicemask, zmask ) 
    1015       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1016          &                  igrdi, igrdj, pseaicen, zseaicel ) 
    1017        
    1018       DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 
    1019           
    1020          iobs = jobs - seaicedatqc%nsurfup 
    1021           
    1022          IF ( kt /= seaicedatqc%mstp(jobs) ) THEN 
    1023              
    1024             IF(lwp) THEN 
    1025                WRITE(numout,*) 
    1026                WRITE(numout,*) ' E R R O R : Observation',              & 
    1027                   &            ' time step is not consistent with the', & 
    1028                   &            ' model time step' 
    1029                WRITE(numout,*) ' =========' 
    1030                WRITE(numout,*) 
    1031                WRITE(numout,*) ' Record  = ', jobs,                & 
    1032                   &            ' kt      = ', kt,                  & 
    1033                   &            ' mstp    = ', seaicedatqc%mstp(jobs), & 
    1034                   &            ' ntyp    = ', seaicedatqc%ntyp(jobs) 
    1035             ENDIF 
    1036             CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' ) 
    1037              
    1038          ENDIF 
    1039           
    1040          zlam = seaicedatqc%rlam(jobs) 
    1041          zphi = seaicedatqc%rphi(jobs) 
    1042           
    1043          ! Get weights to interpolate the model sea ice to the observation point 
    1044          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    1045             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    1046             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    1047           
    1048          ! ... Interpolate the model sea ice to the observation point 
    1049          CALL obs_int_h2d( 1, 1,      & 
    1050             &              zweig, zseaicel(:,:,iobs),  zext ) 
    1051           
    1052          seaicedatqc%rmod(jobs,1) = zext(1) 
    1053           
    1054       END DO 
    1055        
    1056       ! Deallocate the data for interpolation 
    1057       DEALLOCATE( & 
    1058          & igrdi,    & 
    1059          & igrdj,    & 
    1060          & zglam,    & 
    1061          & zgphi,    & 
    1062          & zmask,    & 
    1063          & zseaicel  & 
    1064          & ) 
    1065        
    1066       seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice 
    1067  
    1068    END SUBROUTINE obs_seaice_opt 
    1069  
    1070    SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    1071       &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, & 
    1072       &                    ld_dailyav ) 
    1073       !!----------------------------------------------------------------------- 
    1074       !! 
    1075       !!                     ***  ROUTINE obs_vel_opt  *** 
    1076       !! 
    1077       !! ** Purpose : Compute the model counterpart of velocity profile 
    1078       !!              data by interpolating from the model grid to the  
    1079       !!              observation point. 
    1080       !! 
    1081       !! ** Method  : Linearly interpolate zonal and meridional components of velocity  
    1082       !!              to each observation point using the model values at the corners of  
    1083       !!              the surrounding grid box. The model velocity components are on a  
    1084       !!              staggered C- grid. 
    1085       !! 
    1086       !!    For velocity data from the TAO array, the model equivalent is 
    1087       !!    a daily mean velocity field. So, we first compute 
    1088       !!    the mean, then interpolate only at the end of the day. 
    1089       !! 
    1090       !! ** Action  : 
    1091       !! 
    1092       !! History : 
    1093       !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles 
    1094       !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles 
    1095       !!----------------------------------------------------------------------- 
    1096      
    1097       !! * Modules used 
    1098       USE obs_profiles_def ! Definition of storage space for profile obs. 
    1099  
    1100       IMPLICIT NONE 
    1101  
    1102       !! * Arguments 
    1103       TYPE(obs_prof), INTENT(INOUT) :: & 
    1104          & prodatqc        ! Subset of profile data not failing screening 
    1105       INTEGER, INTENT(IN) :: kt        ! Time step 
    1106       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    1107       INTEGER, INTENT(IN) :: kpj 
    1108       INTEGER, INTENT(IN) :: kpk  
    1109       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    1110                                        !   (kit000-1 = restart time) 
    1111       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header) 
    1112       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    1113       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
    1114       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    1115          & pun,    &    ! Model zonal component of velocity 
    1116          & pvn,    &    ! Model meridional component of velocity 
    1117          & pumask, &    ! Land-sea mask 
    1118          & pvmask       ! Land-sea mask 
    1119       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    1120          & pgdept       ! Model array of depth levels 
    1121       LOGICAL, INTENT(IN) :: ld_dailyav 
    1122           
    1123       !! * Local declarations 
    1124       INTEGER :: ji 
    1125       INTEGER :: jj 
    1126       INTEGER :: jk 
    1127       INTEGER :: jobs 
    1128       INTEGER :: inrc 
    1129       INTEGER :: ipro 
    1130       INTEGER :: idayend 
    1131       INTEGER :: ista 
    1132       INTEGER :: iend 
    1133       INTEGER :: iobs 
    1134       INTEGER, DIMENSION(imaxavtypes) :: & 
    1135          & idailyavtypes 
    1136       REAL(KIND=wp) :: zlam 
    1137       REAL(KIND=wp) :: zphi 
    1138       REAL(KIND=wp) :: zdaystp 
    1139       REAL(KIND=wp), DIMENSION(kpk) :: & 
    1140          & zobsmasku, & 
    1141          & zobsmaskv, & 
    1142          & zobsmask,  & 
    1143          & zobsk,     & 
    1144          & zobs2k 
    1145       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
    1146          & zweigu,zweigv 
    1147       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    1148          & zumask, zvmask, & 
    1149          & zintu, & 
    1150          & zintv, & 
    1151          & zinmu, & 
    1152          & zinmv 
    1153       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    1154          & zglamu, zglamv, & 
    1155          & zgphiu, zgphiv 
    1156       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1157          & igrdiu, & 
    1158          & igrdju, & 
    1159          & igrdiv, & 
    1160          & igrdjv 
    1161  
    1162       !------------------------------------------------------------------------ 
    1163       ! Local initialization  
    1164       !------------------------------------------------------------------------ 
    1165       ! ... Record and data counters 
    1166       inrc = kt - kit000 + 2 
    1167       ipro = prodatqc%npstp(inrc) 
    1168  
    1169       ! Initialize daily mean for first timestep 
    1170       idayend = MOD( kt - kit000 + 1, kdaystp ) 
    1171  
    1172       ! Added kt == 0 test to catch restart case  
    1173       IF ( idayend == 1 .OR. kt == 0) THEN 
    1174          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
    1175          prodatqc%vdmean(:,:,:,1) = 0.0 
    1176          prodatqc%vdmean(:,:,:,2) = 0.0 
    1177       ENDIF 
    1178  
    1179       ! Increment the zonal velocity field for computing daily mean 
    1180       prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:) 
    1181       ! Increment the meridional velocity field for computing daily mean 
    1182       prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:) 
    1183     
    1184       ! Compute the daily mean at the end of day 
    1185       zdaystp = 1.0 / REAL( kdaystp ) 
    1186       IF ( idayend == 0 ) THEN 
    1187          prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp 
    1188          prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp 
    1189       ENDIF 
    1190  
    1191       ! Get the data for interpolation 
    1192       ALLOCATE( & 
    1193          & igrdiu(2,2,ipro),      & 
    1194          & igrdju(2,2,ipro),      & 
    1195          & igrdiv(2,2,ipro),      & 
    1196          & igrdjv(2,2,ipro),      & 
    1197          & zglamu(2,2,ipro), zglamv(2,2,ipro), & 
    1198          & zgphiu(2,2,ipro), zgphiv(2,2,ipro), & 
    1199          & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), & 
    1200          & zintu(2,2,kpk,ipro),  & 
    1201          & zintv(2,2,kpk,ipro)   & 
    1202          & ) 
    1203  
    1204       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    1205          iobs = jobs - prodatqc%nprofup 
    1206          igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    1207          igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    1208          igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    1209          igrdju(1,2,iobs) = prodatqc%mj(jobs,1) 
    1210          igrdiu(2,1,iobs) = prodatqc%mi(jobs,1) 
    1211          igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    1212          igrdiu(2,2,iobs) = prodatqc%mi(jobs,1) 
    1213          igrdju(2,2,iobs) = prodatqc%mj(jobs,1) 
    1214          igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    1215          igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    1216          igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    1217          igrdjv(1,2,iobs) = prodatqc%mj(jobs,2) 
    1218          igrdiv(2,1,iobs) = prodatqc%mi(jobs,2) 
    1219          igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    1220          igrdiv(2,2,iobs) = prodatqc%mi(jobs,2) 
    1221          igrdjv(2,2,iobs) = prodatqc%mj(jobs,2) 
    1222       END DO 
    1223  
    1224       CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu ) 
    1225       CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu ) 
    1226       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask ) 
    1227       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu ) 
    1228  
    1229       CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv ) 
    1230       CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv ) 
    1231       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask ) 
    1232       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv ) 
    1233  
    1234       ! At the end of the day also get interpolated means 
    1235       IF ( idayend == 0 ) THEN 
    1236  
    1237          ALLOCATE( & 
    1238             & zinmu(2,2,kpk,ipro),  & 
    1239             & zinmv(2,2,kpk,ipro)   & 
    1240             & ) 
    1241  
    1242          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, & 
    1243             &                  prodatqc%vdmean(:,:,:,1), zinmu ) 
    1244          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, & 
    1245             &                  prodatqc%vdmean(:,:,:,2), zinmv ) 
    1246  
    1247       ENDIF 
    1248  
    1249 ! loop over observations 
    1250  
    1251       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    1252  
    1253          iobs = jobs - prodatqc%nprofup 
    1254  
    1255          IF ( kt /= prodatqc%mstp(jobs) ) THEN 
    1256              
    1257             IF(lwp) THEN 
    1258                WRITE(numout,*) 
    1259                WRITE(numout,*) ' E R R O R : Observation',              & 
    1260                   &            ' time step is not consistent with the', & 
    1261                   &            ' model time step' 
    1262                WRITE(numout,*) ' =========' 
    1263                WRITE(numout,*) 
    1264                WRITE(numout,*) ' Record  = ', jobs,                    & 
    1265                   &            ' kt      = ', kt,                      & 
    1266                   &            ' mstp    = ', prodatqc%mstp(jobs), & 
    1267                   &            ' ntyp    = ', prodatqc%ntyp(jobs) 
    1268             ENDIF 
    1269             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
    1270          ENDIF 
    1271           
    1272          zlam = prodatqc%rlam(jobs) 
    1273          zphi = prodatqc%rphi(jobs) 
    1274  
    1275          ! Initialize observation masks 
    1276  
    1277          zobsmasku(:) = 0.0 
    1278          zobsmaskv(:) = 0.0 
    1279           
    1280          ! Horizontal weights and vertical mask 
    1281  
    1282          IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    1283  
    1284             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    1285                &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), & 
    1286                &                   zumask(:,:,:,iobs), zweigu, zobsmasku ) 
    1287  
    1288          ENDIF 
    1289  
    1290           
    1291          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    1292  
    1293             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    1294                &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), & 
    1295                &                   zvmask(:,:,:,iobs), zweigv, zobsmasku ) 
    1296  
    1297          ENDIF 
    1298  
    1299          ! Ensure that the vertical mask on u and v are consistent. 
    1300  
    1301          zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) ) 
    1302  
    1303          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    1304  
    1305             zobsk(:) = obfillflt 
    1306  
    1307        IF ( ld_dailyav ) THEN 
    1308  
    1309                IF ( idayend == 0 )  THEN 
    1310                    
    1311                   ! Daily averaged data 
    1312                    
    1313                   CALL obs_int_h2d( kpk, kpk,      & 
    1314                      &              zweigu, zinmu(:,:,:,iobs), zobsk ) 
    1315                    
    1316                    
    1317                ELSE 
    1318                 
    1319                   CALL ctl_stop( ' A nonzero' //     & 
    1320                      &           ' number of U profile data should' // & 
    1321                      &           ' only occur at the end of a given day' ) 
    1322  
    1323                ENDIF 
    1324            
    1325             ELSE  
    1326                 
    1327                ! Point data 
    1328  
    1329                CALL obs_int_h2d( kpk, kpk,      & 
    1330                   &              zweigu, zintu(:,:,:,iobs), zobsk ) 
    1331  
    1332             ENDIF 
    1333  
    1334             !------------------------------------------------------------- 
    1335             ! Compute vertical second-derivative of the interpolating  
    1336             ! polynomial at obs points 
    1337             !------------------------------------------------------------- 
    1338              
    1339             IF ( k1dint == 1 ) THEN 
    1340                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    1341                   &                  pgdept, zobsmask ) 
    1342             ENDIF 
    1343              
    1344             !----------------------------------------------------------------- 
    1345             !  Vertical interpolation to the observation point 
    1346             !----------------------------------------------------------------- 
    1347             ista = prodatqc%npvsta(jobs,1) 
    1348             iend = prodatqc%npvend(jobs,1) 
    1349             CALL obs_int_z1d( kpk,                & 
    1350                & prodatqc%var(1)%mvk(ista:iend),  & 
    1351                & k1dint, iend - ista + 1,         & 
    1352                & prodatqc%var(1)%vdep(ista:iend), & 
    1353                & zobsk, zobs2k,                   & 
    1354                & prodatqc%var(1)%vmod(ista:iend), & 
    1355                & pgdept, zobsmask ) 
    1356  
    1357          ENDIF 
    1358  
    1359          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    1360  
    1361             zobsk(:) = obfillflt 
    1362  
    1363             IF ( ld_dailyav ) THEN 
    1364  
    1365                IF ( idayend == 0 )  THEN 
    1366  
    1367                   ! Daily averaged data 
    1368                    
    1369                   CALL obs_int_h2d( kpk, kpk,      & 
    1370                      &              zweigv, zinmv(:,:,:,iobs), zobsk ) 
    1371                    
    1372                ELSE 
    1373  
    1374                   CALL ctl_stop( ' A nonzero' //     & 
    1375                      &           ' number of V profile data should' // & 
    1376                      &           ' only occur at the end of a given day' ) 
    1377  
    1378                ENDIF 
    1379  
    1380             ELSE 
    1381                 
    1382                ! Point data 
    1383  
    1384                CALL obs_int_h2d( kpk, kpk,      & 
    1385                   &              zweigv, zintv(:,:,:,iobs), zobsk ) 
    1386  
    1387             ENDIF 
    1388  
    1389  
    1390             !------------------------------------------------------------- 
    1391             ! Compute vertical second-derivative of the interpolating  
    1392             ! polynomial at obs points 
    1393             !------------------------------------------------------------- 
    1394              
    1395             IF ( k1dint == 1 ) THEN 
    1396                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    1397                   &                  pgdept, zobsmask ) 
    1398             ENDIF 
    1399              
    1400             !---------------------------------------------------------------- 
    1401             !  Vertical interpolation to the observation point 
    1402             !---------------------------------------------------------------- 
    1403             ista = prodatqc%npvsta(jobs,2) 
    1404             iend = prodatqc%npvend(jobs,2) 
    1405             CALL obs_int_z1d( kpk, & 
    1406                & prodatqc%var(2)%mvk(ista:iend),& 
    1407                & k1dint, iend - ista + 1, & 
    1408                & prodatqc%var(2)%vdep(ista:iend),& 
    1409                & zobsk, zobs2k, & 
    1410                & prodatqc%var(2)%vmod(ista:iend),& 
    1411                & pgdept, zobsmask ) 
    1412  
    1413          ENDIF 
    1414  
    1415       END DO 
    1416   
    1417       ! Deallocate the data for interpolation 
    1418       DEALLOCATE( & 
    1419          & igrdiu, & 
    1420          & igrdju, & 
    1421          & igrdiv, & 
    1422          & igrdjv, & 
    1423          & zglamu, zglamv, & 
    1424          & zgphiu, zgphiv, & 
    1425          & zumask, zvmask, & 
    1426          & zintu, & 
    1427          & zintv  & 
    1428          & ) 
    1429       ! At the end of the day also get interpolated means 
    1430       IF ( idayend == 0 ) THEN 
    1431          DEALLOCATE( & 
    1432             & zinmu,  & 
    1433             & zinmv   & 
    1434             & ) 
    1435       ENDIF 
    1436  
    1437       prodatqc%nprofup = prodatqc%nprofup + ipro  
    1438        
    1439    END SUBROUTINE obs_vel_opt 
     730 
     731      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     732 
     733   END SUBROUTINE obs_surf_opt 
    1440734 
    1441735END MODULE obs_oper 
Note: See TracChangeset for help on using the changeset viewer.