Changeset 12035


Ignore:
Timestamp:
2019-12-03T12:20:53+01:00 (9 months ago)
Author:
kingr
Message:

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

Location:
branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/CONFIG/SHARED/namelist_ref

    r11863 r12035  
    12341234   ln_sstnight = .false.            ! Logical switch for calculating night-time average for SST obs 
    12351235   ln_output_clim = .false.         ! Logical switch for writing climatological values to fdbk files 
     1236   ln_time_mean_sla_bkg = .false.   ! Logical switch for applying time mean of SLA background to remove tidal signal 
    12361237   ln_default_fp_indegs = .true.    ! Logical: T=> averaging footprint is in degrees, F=> in metres 
    12371238   ln_sla_fp_indegs = .true.        ! Logical: T=> averaging footprint is in degrees, F=> in metres 
  • branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r11863 r12035  
    5555   LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
    5656   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology 
     57   LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 
    5758 
    5859   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
     
    6667   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
    6768   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
     69   REAL(wp), PUBLIC :: & 
     70      &        MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 
     71 
    6872 
    6973   INTEGER :: nn_1dint         !: Vertical interpolation method 
     
    238242      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
    239243      LOGICAL :: ltype_clim      ! Local version of ln_output_clim 
     244       
    240245 
    241246      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     
    262267         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    263268         &            ln_sstnight,  ln_output_clim,                   & 
    264          &            ln_default_fp_indegs,                           & 
     269         &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     & 
    265270         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    266271         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     
    427432         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    428433         WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim 
     434         WRITE(numout,*) '             Logical switch for time-mean of SLA        ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg 
    429435      ENDIF 
    430436      !----------------------------------------------------------------------- 
     
    912918               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    913919               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    914                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., & 
     920               &               rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., & 
    915921               &               llnightav(jtype), ltype_clim, clvars ) 
    916922 
     
    10311037      !! * Local declarations 
    10321038      INTEGER :: idaystp           ! Number of timesteps per day 
     1039      INTEGER :: imeanstp          ! Number of timesteps for sla averaging 
    10331040      INTEGER :: jtype             ! Data loop variable 
    10341041      INTEGER :: jvar              ! Variable number 
     
    16241631            ENDIF 
    16251632 
    1626             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1627                &               nit000, idaystp, zsurfvar,               & 
    1628                &               zsurfclim, zsurfmask,                    & 
    1629                &               n2dintsurf(jtype), llnightav(jtype),     & 
    1630                &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    1631                &               lfpindegs(jtype) ) 
     1633            IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND.                 & 
     1634                  &  ln_time_mean_sla_bkg ) THEN 
     1635               !Number of time-steps in meaning period 
     1636               imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 
     1637               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1638                  &               nit000, idaystp, zsurfvar,               & 
     1639                  &               zsurfclim, zsurfmask,                    & 
     1640                  &               n2dintsurf(jtype), llnightav(jtype),     & 
     1641                  &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1642                  &               lfpindegs(jtype), kmeanstp = imeanstp ) 
     1643 
     1644            ELSE 
     1645               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1646                  &               nit000, idaystp, zsurfvar,               & 
     1647                  &               zsurfclim, zsurfmask,                    & 
     1648                  &               n2dintsurf(jtype), llnightav(jtype),     & 
     1649                  &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1650                  &               lfpindegs(jtype) ) 
     1651            ENDIF 
    16321652 
    16331653         END DO 
  • branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r11468 r12035  
    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 surfdataaqc%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  & 
  • branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r11863 r12035  
    3939 
    4040   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, & 
    41       &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
     41      &                     kvars, kextr, kstp, ddobsini, ddobsend, MeanPeriodHours, & 
    4242      &                     ldignmis, ldmod, ldnightav, ldclim, cdvars ) 
    4343      !!--------------------------------------------------------------------- 
     
    7474      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS 
    7575      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS 
     76      REAL(wp), INTENT(IN) :: MeanPeriodHours ! Averaging period in hours 
    7677      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    7778 
     
    280281               inpfiles(jj)%iobsj = -1 
    281282            ENDIF 
     283 
     284            !If the observations are representing a time mean then set the time 
     285            !of the obs to the end of that meaning period relative to the start of the run 
     286            IF ( MeanPeriodHours > 0._wp ) THEN 
     287               DO ji = 1, inpfiles(jj)%nobs 
     288                  inpfiles(jj)%ptim(ji) = & 
     289                        & djulini(jj) + (MeanPeriodHours/24.) 
     290               END DO 
     291            ENDIF    
     292 
    282293            inowin = 0 
    283294            DO ji = 1, inpfiles(jj)%nobs 
Note: See TracChangeset for help on using the changeset viewer.