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 5682 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 – NEMO

Ignore:
Timestamp:
2015-08-12T17:46:45+02:00 (9 years ago)
Author:
mattmartin
Message:

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.

File:
1 edited

Legend:

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

    r5659 r5682  
    1111   !!---------------------------------------------------------------------- 
    1212 
    13    !! * Modules used    
     13   !! * Modules used 
    1414   USE par_kind, ONLY : &         ! Precision variables 
    1515      & wp 
    1616   USE in_out_manager             ! I/O manager 
    1717   USE obs_inter_sup              ! Interpolation support 
    18    USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt 
     18   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt 
    1919      & obs_int_h2d, & 
    2020      & obs_int_h2d_init 
    21    USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt 
     21   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    2222      & obs_int_z1d,    & 
    2323      & obs_int_z1d_spl 
    24    USE obs_const,  ONLY :     & 
    25       & obfillflt      ! Fillvalue    
    26    USE dom_oce,       ONLY : & 
    27       & glamt, glamu, glamv, & 
    28       & gphit, gphiu, gphiv 
    29    USE lib_mpp,       ONLY : & 
     24   USE obs_const,  ONLY :    &    ! Obs fill value 
     25      & obfillflt 
     26   USE dom_oce,       ONLY : &    ! Model lats/lons 
     27      & glamt, & 
     28      & gphit 
     29   USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    3030      & ctl_warn, ctl_stop 
     31   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
     32      & sbc_dcy, nday_qsr 
    3133 
    3234   IMPLICIT NONE 
     
    3537   PRIVATE 
    3638 
    37    PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile observations 
    38       &   obs_surf_opt     ! Compute the model counterpart of surface observations 
    39  
    40    INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 
     39   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
     40      &   obs_surf_opt     ! Compute the model counterpart of surface obs 
     41 
     42   INTEGER, PARAMETER, PUBLIC :: & 
     43      & imaxavtypes = 20   ! Max number of daily avgd obs types 
    4144 
    4245   !!---------------------------------------------------------------------- 
     
    4851CONTAINS 
    4952 
    50    SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    51       &                     pvar1, pvar2, pgdept, pmask1, k1dint, k2dint, & 
    52       &                     kdailyavtypes ) 
     53   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
     54      &                     kit000, kdaystp,                      & 
     55      &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
     56      &                     plam1, plam2, pphi1, pphi2,           & 
     57      &                     k1dint, k2dint, kdailyavtypes ) 
     58 
    5359      !!----------------------------------------------------------------------- 
    5460      !! 
     
    99105      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
    100106      !!----------------------------------------------------------------------- 
    101    
     107 
    102108      !! * Modules used 
    103109      USE obs_profiles_def ! Definition of storage space for profile obs. 
     
    106112 
    107113      !! * Arguments 
    108       TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening 
    109       INTEGER, INTENT(IN) :: kt        ! Time step 
    110       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
     114      TYPE(obs_prof), INTENT(INOUT) :: & 
     115         & prodatqc                  ! Subset of profile data passing QC 
     116      INTEGER, INTENT(IN) :: kt      ! Time step 
     117      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
    111118      INTEGER, INTENT(IN) :: kpj 
    112119      INTEGER, INTENT(IN) :: kpk 
    113       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    114                                        !   (kit000-1 = restart time) 
    115       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header) 
    116       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    117       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
     120      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step 
     121                                     !   (kit000-1 = restart time) 
     122      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header) 
     123      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header) 
     124      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 
    118125      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    119          & pvar1,    &    ! Model field 1 
    120          & pvar2,    &    ! Model field 2 
    121          & pmask1,   &    ! Land-sea mask 1 
    122          & pmask2         ! Land-sea mask 2 
     126         & pvar1,    &               ! Model field 1 
     127         & pvar2,    &               ! Model field 2 
     128         & pmask1,   &               ! Land-sea mask 1 
     129         & pmask2                    ! Land-sea mask 2 
     130      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     131         & plam1,    &               ! Model longitudes for variable 1 
     132         & plam2,    &               ! Model longitudes for variable 2 
     133         & pphi1,    &               ! Model latitudes for variable 1 
     134         & pphi2                     ! Model latitudes for variable 2 
    123135      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    124          & pgdept       ! Model array of depth levels 
     136         & pgdept                    ! Model array of depth levels 
    125137      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    126          & kdailyavtypes! Types for daily averages 
     138         & kdailyavtypes             ! Types for daily averages 
     139 
    127140      !! * Local declarations 
    128141      INTEGER ::   ji 
     
    138151      INTEGER, DIMENSION(imaxavtypes) :: & 
    139152         & idailyavtypes 
     153      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     154         & igrdi1, & 
     155         & igrdi2, & 
     156         & igrdj1, & 
     157         & igrdj2 
    140158      REAL(KIND=wp) :: zlam 
    141159      REAL(KIND=wp) :: zphi 
     
    144162         & zobsmask1, & 
    145163         & zobsmask2, & 
    146          & zobsmask, & 
    147164         & zobsk,    & 
    148165         & zobs2k 
     
    162179         & zgphi1, & 
    163180         & zgphi2 
    164       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    165          & igrdi1, & 
    166          & igrdi2, & 
    167          & igrdj1, & 
    168          & igrdj2 
     181      LOGICAL :: ld_dailyav 
    169182 
    170183      !------------------------------------------------------------------------ 
    171184      ! Local initialization  
    172185      !------------------------------------------------------------------------ 
    173       ! ... Record and data counters 
     186      ! Record and data counters 
    174187      inrc = kt - kit000 + 2 
    175188      ipro = prodatqc%npstp(inrc) 
    176   
     189 
    177190      ! Daily average types 
     191      ld_dailyav = .FALSE. 
    178192      IF ( PRESENT(kdailyavtypes) ) THEN 
    179193         idailyavtypes(:) = kdailyavtypes(:) 
     194         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 
    180195      ELSE 
    181196         idailyavtypes(:) = -1 
    182197      ENDIF 
    183198 
    184       ! Initialize daily mean for first timestep 
     199      ! Daily means are calculated for values over timesteps: 
     200      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 
    185201      idayend = MOD( kt - kit000 + 1, kdaystp ) 
    186202 
    187       ! Added kt == 0 test to catch restart case  
    188       IF ( idayend == 1 .OR. kt == 0) THEN 
    189          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
     203      IF ( ld_dailyav ) THEN 
     204 
     205         ! Initialize daily mean for first timestep of the day 
     206         IF ( idayend == 1 .OR. kt == 0 ) THEN 
     207            DO jk = 1, jpk 
     208               DO jj = 1, jpj 
     209                  DO ji = 1, jpi 
     210                     prodatqc%vdmean(ji,jj,jk,1) = 0.0 
     211                     prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     212                  END DO 
     213               END DO 
     214            END DO 
     215         ENDIF 
     216 
    190217         DO jk = 1, jpk 
    191218            DO jj = 1, jpj 
    192219               DO ji = 1, jpi 
    193                   prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    194                   prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     220                  ! Increment field 1 for computing daily mean 
     221                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
     222                     &                        + pvar1(ji,jj,jk) 
     223                  ! Increment field 2 for computing daily mean 
     224                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     225                     &                        + pvar2(ji,jj,jk) 
    195226               END DO 
    196227            END DO 
    197228         END DO 
    198       ENDIF 
    199  
    200       DO jk = 1, jpk 
    201          DO jj = 1, jpj 
    202             DO ji = 1, jpi 
    203                ! Increment field 1 for computing daily mean 
    204                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    205                   &                        + pvar1(ji,jj,jk) 
    206                ! Increment field 2 for computing daily mean 
    207                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    208                   &                        + pvar2(ji,jj,jk) 
    209             END DO 
    210          END DO 
    211       END DO 
    212     
    213       ! Compute the daily mean at the end of day 
    214       zdaystp = 1.0 / REAL( kdaystp ) 
    215       IF ( idayend == 0 ) THEN 
    216          DO jk = 1, jpk 
    217             DO jj = 1, jpj 
    218                DO ji = 1, jpi 
    219                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    220                      &                        * zdaystp 
    221                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    222                   &                           * zdaystp 
     229 
     230         ! Compute the daily mean at the end of day 
     231         zdaystp = 1.0 / REAL( kdaystp ) 
     232         IF ( idayend == 0 ) THEN 
     233            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 
     234            CALL FLUSH(numout) 
     235            DO jk = 1, jpk 
     236               DO jj = 1, jpj 
     237                  DO ji = 1, jpi 
     238                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
     239                        &                        * zdaystp 
     240                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     241                        &                        * zdaystp 
     242                  END DO 
    223243               END DO 
    224244            END DO 
    225          END DO 
     245         ENDIF 
     246 
    226247      ENDIF 
    227248 
     
    262283      END DO 
    263284 
    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 ) 
     285      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, plam1, zglam1 ) 
     286      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, pphi1, zgphi1 ) 
    266287      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 
    267288      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
    268289       
    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 ) 
     290      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, plam2, zglam2 ) 
     291      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, pphi2, zgphi2 ) 
     292      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
    272293      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    273294 
    274295      ! At the end of the day also get interpolated means 
    275       IF ( idayend == 0 ) THEN 
     296      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    276297 
    277298         ALLOCATE( & 
     
    292313 
    293314         IF ( kt /= prodatqc%mstp(jobs) ) THEN 
    294              
     315 
    295316            IF(lwp) THEN 
    296317               WRITE(numout,*) 
     
    307328            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
    308329         ENDIF 
    309           
     330 
    310331         zlam = prodatqc%rlam(jobs) 
    311332         zphi = prodatqc%rphi(jobs) 
    312           
     333 
    313334         ! Horizontal weights and vertical mask 
    314335 
     
    333354            zobsk(:) = obfillflt 
    334355 
    335        IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
     356            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    336357 
    337358               IF ( idayend == 0 )  THEN 
    338                    
    339359                  ! Daily averaged data 
    340                    
    341360                  CALL obs_int_h2d( kpk, kpk,      & 
    342361                     &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    343                    
    344                    
    345                ELSE 
    346                 
    347                   CALL ctl_stop( ' A nonzero' //     & 
    348                      &           ' number of average profile data should' // & 
    349                      &           ' only occur at the end of a given day' ) 
    350362 
    351363               ENDIF 
    352            
     364 
    353365            ELSE  
    354                 
     366 
    355367               ! Point data 
    356  
    357368               CALL obs_int_h2d( kpk, kpk,      & 
    358369                  &              zweig1, zint1(:,:,:,iobs), zobsk ) 
     
    364375            ! polynomial at obs points 
    365376            !------------------------------------------------------------- 
    366              
     377 
    367378            IF ( k1dint == 1 ) THEN 
    368379               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    369                   &                  pgdept, zobsmask ) 
     380                  &                  pgdept, zobsmask1 ) 
    370381            ENDIF 
    371              
     382 
    372383            !----------------------------------------------------------------- 
    373384            !  Vertical interpolation to the observation point 
     
    381392               & zobsk, zobs2k,                   & 
    382393               & prodatqc%var(1)%vmod(ista:iend), & 
    383                & pgdept, zobsmask ) 
     394               & pgdept, zobsmask1 ) 
    384395 
    385396         ENDIF 
     
    394405 
    395406                  ! Daily averaged data 
    396                    
    397407                  CALL obs_int_h2d( kpk, kpk,      & 
    398408                     &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    399                    
    400                ELSE 
    401  
    402                   CALL ctl_stop( ' A nonzero' //     & 
    403                      &           ' number of average profile data should' // & 
    404                      &           ' only occur at the end of a given day' ) 
    405409 
    406410               ENDIF 
    407411 
    408412            ELSE 
    409                 
     413 
    410414               ! Point data 
    411  
    412415               CALL obs_int_h2d( kpk, kpk,      & 
    413416                  &              zweig2, zint2(:,:,:,iobs), zobsk ) 
     
    420423            ! polynomial at obs points 
    421424            !------------------------------------------------------------- 
    422              
     425 
    423426            IF ( k1dint == 1 ) THEN 
    424427               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    425                   &                  pgdept, zobsmask ) 
     428                  &                  pgdept, zobsmask2 ) 
    426429            ENDIF 
    427              
     430 
    428431            !---------------------------------------------------------------- 
    429432            !  Vertical interpolation to the observation point 
     
    437440               & zobsk, zobs2k, & 
    438441               & prodatqc%var(2)%vmod(ista:iend),& 
    439                & pgdept, zobsmask ) 
     442               & pgdept, zobsmask2 ) 
    440443 
    441444         ENDIF 
    442445 
    443446      END DO 
    444   
     447 
    445448      ! Deallocate the data for interpolation 
    446449      DEALLOCATE( & 
     
    458461         & zint2   & 
    459462         & ) 
     463 
    460464      ! At the end of the day also get interpolated means 
    461       IF ( idayend == 0 ) THEN 
     465      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    462466         DEALLOCATE( & 
    463467            & zinm1,  & 
     
    467471 
    468472      prodatqc%nprofup = prodatqc%nprofup + ipro  
    469        
     473 
    470474   END SUBROUTINE obs_prof_opt 
    471475 
    472    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, kit000, & 
    473       &                    psurf, psurfmask, k2dint, ld_nightav ) 
     476   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         & 
     477      &                    kit000, kdaystp, psurf, psurfmask, & 
     478      &                    k2dint, ldnightav ) 
     479 
    474480      !!----------------------------------------------------------------------- 
    475481      !! 
     
    491497      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    492498      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    493       !!   
     499      !! 
    494500      !! 
    495501      !! ** Action  : 
     
    499505      !!      ! 15-02 (M. Martin) Combined routine for surface types 
    500506      !!----------------------------------------------------------------------- 
    501    
     507 
    502508      !! * Modules used 
    503509      USE obs_surf_def  ! Definition of storage space for surface observations 
     
    506512 
    507513      !! * Arguments 
    508       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc     ! Subset of surface data not failing screening 
    509       INTEGER, INTENT(IN) :: kt      ! Time step 
    510       INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
     514      TYPE(obs_surf), INTENT(INOUT) :: & 
     515         & surfdataqc                  ! Subset of surface data passing QC 
     516      INTEGER, INTENT(IN) :: kt        ! Time step 
     517      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    511518      INTEGER, INTENT(IN) :: kpj 
    512       INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
    513                                       !   (kit000-1 = restart time) 
    514       INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    515       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    516          & psurf,  &    ! Model surface field 
    517          & psurfmask    ! Land-sea mask 
    518           
     519      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
     520                                       !   (kit000-1 = restart time) 
     521      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
     522      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     523      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     524         & psurf,  &                   ! Model surface field 
     525         & psurfmask                   ! Land-sea mask 
     526      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     527 
    519528      !! * Local declarations 
    520529      INTEGER :: ji 
     
    524533      INTEGER :: isurf 
    525534      INTEGER :: iobs 
    526       REAL(KIND=wp) :: zlam 
    527       REAL(KIND=wp) :: zphi 
    528       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    529       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    530          & zweig 
    531       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    532          & zmask, & 
    533          & zsurf, & 
    534          & zglam, & 
    535          & zgphi 
     535      INTEGER :: idayend 
    536536      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    537537         & igrdi, & 
    538538         & igrdj 
     539      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     540         & icount_night,      & 
     541         & imask_night 
     542      REAL(wp) :: zlam 
     543      REAL(wp) :: zphi 
     544      REAL(wp), DIMENSION(1) :: zext, zobsmask 
     545      REAL(wp) :: zdaystp 
     546      REAL(wp), DIMENSION(2,2,1) :: & 
     547         & zweig 
     548      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     549         & zmask,  & 
     550         & zsurf,  & 
     551         & zsurfm, & 
     552         & zglam,  & 
     553         & zgphi 
     554      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     555         & zintmp,  & 
     556         & zouttmp, & 
     557         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    539558 
    540559      !------------------------------------------------------------------------ 
    541560      ! Local initialization  
    542561      !------------------------------------------------------------------------ 
    543       ! ... Record and data counters 
     562      ! Record and data counters 
    544563      inrc = kt - kit000 + 2 
    545564      isurf = surfdataqc%nsstp(inrc) 
    546        
    547       IF ( ld_nightav ) THEN 
     565 
     566      IF ( ldnightav ) THEN 
    548567 
    549568      ! Initialize array for night mean 
    550          IF ( kt .EQ. 0 ) THEN 
     569         IF ( kt == 0 ) THEN 
    551570            ALLOCATE ( icount_night(kpi,kpj) ) 
    552571            ALLOCATE ( imask_night(kpi,kpj) ) 
     
    557576         ENDIF 
    558577 
    559          ! Initialize daily mean for first timestep 
     578         ! Night-time means are calculated for night-time values over timesteps: 
     579         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 
    560580         idayend = MOD( kt - kit000 + 1, kdaystp ) 
    561581 
    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 
     582         ! Initialize night-time mean for first timestep of the day 
     583         IF ( idayend == 1 .OR. kt == 0 ) THEN 
    565584            DO jj = 1, jpj 
    566585               DO ji = 1, jpi 
     
    580599               ! Increment the temperature field for computing night mean and counter 
    581600               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) 
     601                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 
     602               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj) 
     603               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj) 
    585604            END DO 
    586605         END DO 
    587     
    588          ! Compute the daily mean at the end of day 
    589  
     606 
     607         ! Compute the night-time mean at the end of the day 
    590608         zdaystp = 1.0 / REAL( kdaystp ) 
    591  
    592          IF ( idayend == 0 ) THEN  
     609         IF ( idayend == 0 ) THEN 
     610            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 
    593611            DO jj = 1, jpj 
    594612               DO ji = 1, jpi 
    595613                  ! Test if "no night" point 
    596                   IF ( icount_night(ji,jj) .NE. 0 ) THEN 
     614                  IF ( icount_night(ji,jj) > 0 ) THEN 
    597615                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
    598                        &                        / icount_night(ji,jj)  
     616                       &                        / REAL( icount_night(ji,jj) ) 
    599617                  ELSE 
     618                     !At locations where there is no night (e.g. poles), 
     619                     ! calculate daily mean instead of night-time mean. 
    600620                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
    601621                  ENDIF 
     
    616636         & zsurf(2,2,isurf)  & 
    617637         & ) 
    618        
     638 
    619639      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    620640         iobs = jobs - surfdataqc%nsurfup 
     
    639659 
    640660      ! At the end of the day get interpolated means 
    641       IF ( idayend == 0 .AND. ld_nightav ) THEN 
     661      IF ( idayend == 0 .AND. ldnightav ) THEN 
    642662 
    643663         ALLOCATE( & 
     
    656676 
    657677         IF ( kt /= surfdataqc%mstp(jobs) ) THEN 
    658              
     678 
    659679            IF(lwp) THEN 
    660680               WRITE(numout,*) 
     
    669689                  &            ' ntyp    = ', surfdataqc%ntyp(jobs) 
    670690            ENDIF 
    671             CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 
    672              
    673          ENDIF 
    674           
     691            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 
     692 
     693         ENDIF 
     694 
    675695         zlam = surfdataqc%rlam(jobs) 
    676696         zphi = surfdataqc%rphi(jobs) 
     
    680700            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    681701            &                   zmask(:,:,iobs), zweig, zobsmask ) 
    682           
    683702 
    684703         ! 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  
     704         IF ( ldnightav .AND. idayend == 0 ) THEN 
     705            ! Night-time averaged data 
     706            CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
    697707         ELSE 
    698  
    699             CALL obs_int_h2d( 1, 1,      & 
    700                &              zweig, zsurf(:,:,iobs),  zext ) 
    701  
    702          ENDIF 
    703           
    704          IF ( surfdataqc%nextr == 2 ) THEN 
     708            CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     709         ENDIF 
     710 
     711         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
    705712            ! ... Remove the MDT from the SSH at the observation point to get the SLA 
    706713            surfdataqc%rext(jobs,1) = zext(1) 
     
    722729         & ) 
    723730 
    724       ! At the end of the day also get interpolated means 
    725       IF ( idayend == 0 .AND. ld_nightav ) THEN 
     731      ! At the end of the day also deallocate night-time mean array 
     732      IF ( idayend == 0 .AND. ldnightav ) THEN 
    726733         DEALLOCATE( & 
    727734            & zsurfm  & 
     
    734741 
    735742END MODULE obs_oper 
    736  
Note: See TracChangeset for help on using the changeset viewer.