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 12140 for branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 – NEMO

Ignore:
Timestamp:
2019-12-10T12:38:10+01:00 (4 years ago)
Author:
kingr
Message:

Added option to remove tides from SLA bkg by taking average over 24h50m.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r11468 r12140  
    542542      &                     kit000, kdaystp, psurf, pclim, psurfmask,   & 
    543543      &                     k2dint, ldnightav, plamscl, pphiscl, & 
    544       &                     lindegrees ) 
     544      &                     lindegrees, kmeanstp ) 
    545545 
    546546      !!----------------------------------------------------------------------- 
     
    592592      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
    593593      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     594      INTEGER, INTENT(IN), OPTIONAL :: & 
     595                             kmeanstp  ! Number of time steps for the time meaning 
     596                                       ! Averaging is triggered if present and greater than one                     
    594597      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    595598         & psurf,  &                   ! Model surface field 
     
    613616      INTEGER :: imodi, imodj 
    614617      INTEGER :: idayend 
     618      INTEGER :: imeanend 
    615619      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    616620         & igrdi,   & 
     
    625629      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 
    626630      REAL(wp) :: zdaystp 
     631      REAL(wp) :: zmeanstp 
    627632      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    628633         & zweig,  & 
     
    641646         & zouttmp, & 
    642647         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
     648 
     649      LOGICAL :: l_timemean 
    643650          
    644651      !------------------------------------------------------------------------ 
     
    649656      isurf = surfdataqc%nsstp(inrc) 
    650657 
     658      l_timemean = .FALSE. 
     659      IF ( PRESENT( kmeanstp ) ) THEN 
     660         IF ( kmeanstp > 1 ) l_timemean = .TRUE. 
     661      ENDIF 
     662 
    651663      ! Work out the maximum footprint size for the  
    652664      ! interpolation/averaging in model grid-points - has to be even. 
    653665 
    654666      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     667 
     668 
     669      IF ( l_timemean ) THEN 
     670         ! Initialize time mean for first timestep 
     671         imeanend = MOD( kt - kit000 + 1, kmeanstp ) 
     672         IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend 
     673 
     674         ! Added kt == 0 test to catch restart case  
     675         IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN 
     676            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 
     677            DO jj = 1, jpj 
     678               DO ji = 1, jpi 
     679                  surfdataqc%vdmean(ji,jj) = 0.0 
     680               END DO 
     681            END DO 
     682         ENDIF 
     683 
     684         ! On each time-step, increment the field for computing time mean 
     685         IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt 
     686         DO jj = 1, jpj 
     687            DO ji = 1, jpi 
     688               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     689                  &                        + psurf(ji,jj) 
     690            END DO 
     691         END DO 
     692 
     693         ! Compute the time mean at the end of time period 
     694         IF ( imeanend == 0 ) THEN 
     695            zmeanstp = 1.0 / REAL( kmeanstp ) 
     696            IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp 
     697            DO jj = 1, jpj 
     698               DO ji = 1, jpi 
     699                  surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     700                     &                       * zmeanstp 
     701               END DO 
     702            END DO 
     703         ENDIF 
     704      ENDIF !l_timemean 
    655705 
    656706 
     
    770820      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    771821         &                  igrdi, igrdj, psurfmask, zmask ) 
    772       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    773          &                  igrdi, igrdj, psurf, zsurf ) 
     822 
     823      ! At the end of the averaging period get interpolated means 
     824      IF ( l_timemean ) THEN 
     825         IF ( imeanend == 0 ) THEN 
     826            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) ) 
     827            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 
     828            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     829               &                  igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm ) 
     830         ENDIF 
     831      ELSE 
     832         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     833            &                  igrdi, igrdj, psurf, zsurf ) 
     834      ENDIF 
    774835 
    775836      IF ( k2dint > 4 ) THEN          
     
    823884         zphi = surfdataqc%rphi(jobs) 
    824885 
    825          IF ( ldnightav .AND. idayend == 0 ) THEN 
    826             ! Night-time averaged data 
     886         IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN 
     887            ! Night-time or N=kmeanstp timestep averaged data 
    827888            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    828889         ELSE 
     
    830891         ENDIF 
    831892 
    832          IF ( k2dint <= 4 ) THEN 
    833  
    834             ! Get weights to interpolate the model value to the observation point 
    835             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    836                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    837                &                   zmask(:,:,iobs), zweig, zobsmask ) 
    838  
    839             ! Interpolate the model value to the observation point  
    840             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
    841  
    842             IF ( surfdataqc%lclim ) THEN   
    843                CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
     893         IF (   ( .NOT. l_timemean ) .OR. &  
     894             &  ( l_timemean .AND. imeanend == 0) ) THEN 
     895            IF ( k2dint <= 4 ) THEN 
     896 
     897               ! Get weights to interpolate the model value to the observation point 
     898               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     899                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     900                  &                   zmask(:,:,iobs), zweig, zobsmask ) 
     901 
     902               ! Interpolate the model value to the observation point  
     903               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     904 
     905               IF ( surfdataqc%lclim ) THEN   
     906                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
     907               ENDIF 
     908 
     909 
     910            ELSE 
     911 
     912               ! Get weights to average the model SLA to the observation footprint 
     913               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     914                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     915                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     916                  &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     917                  &                   lindegrees, zweig, zobsmask ) 
     918 
     919               ! Average the model SST to the observation footprint 
     920               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     921                  &              zweig, zsurftmp(:,:,iobs),  zext ) 
     922 
     923               IF ( surfdataqc%lclim ) THEN   
     924                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     925                     &              zweig, zclim(:,:,iobs),  zclm ) 
     926               ENDIF 
     927 
    844928            ENDIF 
    845929 
    846  
    847          ELSE 
    848  
    849             ! Get weights to average the model SLA to the observation footprint 
    850             CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
    851                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    852                &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
    853                &                   zmask(:,:,iobs), plamscl, pphiscl, & 
    854                &                   lindegrees, zweig, zobsmask ) 
    855  
    856             ! Average the model SST to the observation footprint 
    857             CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
    858                &              zweig, zsurftmp(:,:,iobs),  zext ) 
    859  
    860             IF ( surfdataqc%lclim ) THEN   
    861                CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
    862                   &              zweig, zclim(:,:,iobs),  zclm ) 
     930            IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
     931               ! ... Remove the MDT from the SSH at the observation point to get the SLA 
     932               surfdataqc%rext(jobs,1) = zext(1) 
     933               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
     934            ELSE 
     935               surfdataqc%rmod(jobs,1) = zext(1) 
    863936            ENDIF 
    864937 
    865          ENDIF 
    866  
    867          IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
    868             ! ... Remove the MDT from the SSH at the observation point to get the SLA 
    869             surfdataqc%rext(jobs,1) = zext(1) 
    870             surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
    871          ELSE 
    872             surfdataqc%rmod(jobs,1) = zext(1) 
    873          ENDIF 
    874           
    875          IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
    876           
    877          IF ( zext(1) == obfillflt ) THEN 
    878             ! If the observation value is a fill value, set QC flag to bad 
    879             surfdataqc%nqc(jobs) = 4 
     938            IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
     939 
     940            IF ( zext(1) == obfillflt ) THEN 
     941               ! If the observation value is a fill value, set QC flag to bad 
     942               surfdataqc%nqc(jobs) = 4 
     943            ENDIF          
    880944         ENDIF 
    881945 
     
    901965 
    902966      ! At the end of the day also deallocate night-time mean array 
    903       IF ( idayend == 0 .AND. ldnightav ) THEN 
     967      IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN 
    904968         DEALLOCATE( & 
    905969            & zsurfm  & 
Note: See TracChangeset for help on using the changeset viewer.