Changeset 12380


Ignore:
Timestamp:
2020-02-13T15:10:57+01:00 (12 days ago)
Author:
dcarneir
Message:

Updating obs_oper with obs_oper version in the apps trunk

Location:
branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM
Files:
5 edited

Legend:

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

    r11929 r12380  
    12351235   ln_sstnight = .false.            ! Logical switch for calculating night-time average for SST obs 
    12361236   ln_output_clim = .false.         ! Logical switch for writing climatological values to fdbk files 
     1237   ln_time_mean_sla_bkg = .false.   ! Logical switch for applying time mean of SLA background to remove tidal signal 
    12371238   ln_default_fp_indegs = .true.    ! Logical: T=> averaging footprint is in degrees, F=> in metres 
    12381239   ln_sla_fp_indegs = .true.        ! Logical: T=> averaging footprint is in degrees, F=> in metres 
  • branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r9188 r12380  
    4545   USE agrif_opa_interp ! agrif 
    4646#endif 
    47 #if defined key_asminc    
    48    USE asminc          ! Assimilation increment 
    49 #endif 
    5047 
    5148   IMPLICIT NONE 
     
    457454      ELSE 
    458455         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   & 
    459                 &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       ) 
    460       ENDIF 
    461 #if defined key_asminc 
    462       !                                         ! Include the IAU weighted SSH increment 
    463       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    464          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    465       ENDIF 
    466 #endif 
    467       !                                   !* Fill boundary data arrays for AGRIF 
    468       !                                   ! ------------------------------------ 
     456                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
     457      ENDIF 
     458      !                                   !* Fill boundary data arrays with AGRIF 
     459      !                                   ! ------------------------------------- 
    469460#if defined key_agrif 
    470461         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) 
  • branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r12364 r12380  
    5757   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
    5858   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology 
     59   LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 
    5960 
    6061   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
     
    7071   REAL(wp) :: rn_sit_avglamscl     !: E/W diameter of SIT observation footprint 
    7172   REAL(wp) :: rn_sit_avgphiscl     !: N/S diameter of SIT observation footprint 
     73   REAL(wp), PUBLIC :: & 
     74      &        MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 
     75 
    7276 
    7377   INTEGER :: nn_1dint         !: Vertical interpolation method 
     
    274278         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    275279         &            ln_sstnight,  ln_output_clim,                   & 
    276          &            ln_default_fp_indegs,                           & 
     280         &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     & 
    277281         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    278282         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     
    448452         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    449453         WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim 
     454         WRITE(numout,*) '             Logical switch for time-mean of SLA        ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg 
    450455      ENDIF 
    451456      !----------------------------------------------------------------------- 
     
    952957               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    953958               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    954                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., & 
    955                &               llnightav(jtype), ltype_clim, clvars ) 
     959               &               rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., & 
     960               &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 
    956961 
    957962            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject, ln_seaicetypes ) 
     
    10771082      !! * Local declarations 
    10781083      INTEGER :: idaystp           ! Number of timesteps per day 
     1084      INTEGER :: imeanstp          ! Number of timesteps for sla averaging 
    10791085      INTEGER :: jtype             ! Data loop variable 
    10801086      INTEGER :: jvar              ! Variable number 
     
    16861692            ENDIF 
    16871693 
    1688             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1689                &               nit000, idaystp, zsurfvar,               & 
    1690                &               zsurfclim, zsurfmask,                    & 
    1691                &               n2dintsurf(jtype), llnightav(jtype),     & 
    1692                &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    1693                &               lfpindegs(jtype) ) 
    1694  
    1695  
    1696            ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice 
    1697            ! Thickness") 
     1694            IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND.                 & 
     1695                  &  ln_time_mean_sla_bkg ) THEN 
     1696               !Number of time-steps in meaning period 
     1697               imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 
     1698               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1699                  &               nit000, idaystp, zsurfvar,               & 
     1700                  &               zsurfclim, zsurfmask,                    & 
     1701                  &               n2dintsurf(jtype), llnightav(jtype),     & 
     1702                  &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1703                  &               lfpindegs(jtype), kmeanstp = imeanstp ) 
     1704 
     1705 
     1706            ELSE 
     1707               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1708                  &               nit000, idaystp, zsurfvar,               & 
     1709                  &               zsurfclim, zsurfmask,                    & 
     1710                  &               n2dintsurf(jtype), llnightav(jtype),     & 
     1711                  &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1712                  &               lfpindegs(jtype) ) 
     1713            ENDIF 
     1714 
     1715            ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice 
     1716            ! Thickness") 
    16981717            IF ( TRIM(surfdataqc(jtype)%cvars(1)) == 'FBD' ) THEN 
    16991718                 surfdata(jtype)%cvars(1) = 'SIT' 
  • branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r11932 r12380  
    546546      &                     kit000, kdaystp, psurf, pclim, psurfmask,   & 
    547547      &                     k2dint, ldnightav, plamscl, pphiscl, & 
    548       &                     lindegrees ) 
     548      &                     lindegrees, kmeanstp ) 
    549549 
    550550      !!----------------------------------------------------------------------- 
     
    596596      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
    597597      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     598      INTEGER, INTENT(IN), OPTIONAL :: & 
     599                             kmeanstp  ! Number of time steps for the time meaning 
     600                                       ! Averaging is triggered if present and greater than one                     
    598601      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    599602         & psurf,  &                   ! Model surface field 
     
    617620      INTEGER :: imodi, imodj 
    618621      INTEGER :: idayend 
     622      INTEGER :: imeanend 
    619623      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    620624         & igrdi,   & 
     
    629633      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 
    630634      REAL(wp) :: zdaystp 
     635      REAL(wp) :: zmeanstp 
    631636      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    632637         & zweig,  & 
     
    645650         & zouttmp, & 
    646651         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
     652 
     653      LOGICAL :: l_timemean 
    647654          
    648655      !------------------------------------------------------------------------ 
     
    653660      isurf = surfdataqc%nsstp(inrc) 
    654661 
     662      l_timemean = .FALSE. 
     663      IF ( PRESENT( kmeanstp ) ) THEN 
     664         IF ( kmeanstp > 1 ) l_timemean = .TRUE. 
     665      ENDIF 
     666 
    655667      ! Work out the maximum footprint size for the  
    656668      ! interpolation/averaging in model grid-points - has to be even. 
    657669 
    658670      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     671 
     672 
     673      IF ( l_timemean ) THEN 
     674         ! Initialize time mean for first timestep 
     675         imeanend = MOD( kt - kit000 + 1, kmeanstp ) 
     676         IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend 
     677 
     678         ! Added kt == 0 test to catch restart case  
     679         IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN 
     680            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 
     681            DO jj = 1, jpj 
     682               DO ji = 1, jpi 
     683                  surfdataqc%vdmean(ji,jj) = 0.0 
     684               END DO 
     685            END DO 
     686         ENDIF 
     687 
     688         ! On each time-step, increment the field for computing time mean 
     689         IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt 
     690         DO jj = 1, jpj 
     691            DO ji = 1, jpi 
     692               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     693                  &                        + psurf(ji,jj) 
     694            END DO 
     695         END DO 
     696 
     697         ! Compute the time mean at the end of time period 
     698         IF ( imeanend == 0 ) THEN 
     699            zmeanstp = 1.0 / REAL( kmeanstp ) 
     700            IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp 
     701            DO jj = 1, jpj 
     702               DO ji = 1, jpi 
     703                  surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     704                     &                       * zmeanstp 
     705               END DO 
     706            END DO 
     707         ENDIF 
     708      ENDIF !l_timemean 
    659709 
    660710 
     
    774824      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    775825         &                  igrdi, igrdj, psurfmask, zmask ) 
    776       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    777          &                  igrdi, igrdj, psurf, zsurf ) 
     826 
     827      ! At the end of the averaging period get interpolated means 
     828      IF ( l_timemean ) THEN 
     829         IF ( imeanend == 0 ) THEN 
     830            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) ) 
     831            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 
     832            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     833               &                  igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm ) 
     834         ENDIF 
     835      ELSE 
     836         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     837            &                  igrdi, igrdj, psurf, zsurf ) 
     838      ENDIF 
    778839 
    779840      IF ( k2dint > 4 ) THEN          
     
    827888         zphi = surfdataqc%rphi(jobs) 
    828889 
    829          IF ( ldnightav .AND. idayend == 0 ) THEN 
    830             ! Night-time averaged data 
     890         IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN 
     891            ! Night-time or N=kmeanstp timestep averaged data 
    831892            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    832893         ELSE 
     
    834895         ENDIF 
    835896 
    836          IF ( k2dint <= 4 ) THEN 
    837  
    838             ! Get weights to interpolate the model value 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 value to the observation point  
    844             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
    845  
    846             IF ( surfdataqc%lclim ) THEN   
    847                CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
     897         IF (   ( .NOT. l_timemean ) .OR. &  
     898             &  ( l_timemean .AND. imeanend == 0) ) THEN 
     899            IF ( k2dint <= 4 ) THEN 
     900 
     901               ! Get weights to interpolate the model value to the observation point 
     902               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     903                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     904                  &                   zmask(:,:,iobs), zweig, zobsmask ) 
     905 
     906               ! Interpolate the model value to the observation point  
     907               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     908 
     909               IF ( surfdataqc%lclim ) THEN   
     910                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
     911               ENDIF 
     912 
     913 
     914            ELSE 
     915 
     916               ! Get weights to average the model SLA to the observation footprint 
     917               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     918                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     919                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     920                  &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     921                  &                   lindegrees, zweig, zobsmask ) 
     922 
     923               ! Average the model SST to the observation footprint 
     924               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     925                  &              zweig, zsurftmp(:,:,iobs),  zext ) 
     926 
     927               IF ( surfdataqc%lclim ) THEN   
     928                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     929                     &              zweig, zclim(:,:,iobs),  zclm ) 
     930               ENDIF 
     931 
    848932            ENDIF 
    849933 
    850  
    851          ELSE 
    852  
    853             ! Get weights to average the model SLA to the observation footprint 
    854             CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
    855                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    856                &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
    857                &                   zmask(:,:,iobs), plamscl, pphiscl, & 
    858                &                   lindegrees, zweig, zobsmask ) 
    859  
    860             ! Average the model SST to the observation footprint 
    861             CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
    862                &              zweig, zsurftmp(:,:,iobs),  zext ) 
    863  
    864             IF ( surfdataqc%lclim ) THEN   
    865                CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
    866                   &              zweig, zclim(:,:,iobs),  zclm ) 
     934            IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
     935               ! ... Remove the MDT from the SSH at the observation point to get the SLA 
     936               surfdataqc%rext(jobs,1) = zext(1) 
     937               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
     938            ELSE 
     939               surfdataqc%rmod(jobs,1) = zext(1) 
    867940            ENDIF 
    868941 
    869          ENDIF 
    870  
    871          IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
    872             ! ... Remove the MDT from the SSH at the observation point to get the SLA 
    873             surfdataqc%rext(jobs,1) = zext(1) 
    874             surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
    875          ELSE 
    876             surfdataqc%rmod(jobs,1) = zext(1) 
     942            IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
     943 
     944            IF ( zext(1) == obfillflt ) THEN 
     945               ! If the observation value is a fill value, set QC flag to bad 
     946               surfdataqc%nqc(jobs) = 4 
     947            ENDIF          
    877948         ENDIF 
    878949          
     
    900971            surfdataqc%nqc(jobs) = 4 
    901972         ENDIF 
    902  
     973          
    903974      END DO 
    904975 
     
    922993 
    923994      ! At the end of the day also deallocate night-time mean array 
    924       IF ( idayend == 0 .AND. ldnightav ) THEN 
     995      IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN 
    925996         DEALLOCATE( & 
    926997            & zsurfm  & 
  • branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r11932 r12380  
    3939 
    4040   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, & 
    41       &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    42       &                     ldignmis, ldmod, ldnightav, ldclim, cdvars ) 
     41      &                     kvars, kextr, kstp, ddobsini, ddobsend, MeanPeriodHours, & 
     42      &                     ldignmis, ldmod, ldnightav, ldclim, ln_time_mean_sla_bkg, cdvars ) 
    4343      !!--------------------------------------------------------------------- 
    4444      !! 
     
    7272      LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average 
    7373      LOGICAL, INTENT(IN) :: ldclim     ! Will include climatology at obs points. 
     74      LOGICAL, INTENT(IN) :: ln_time_mean_sla_bkg     ! Will reset times to end of averaging period. 
    7475      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS 
    7576      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS 
     77      REAL(wp), INTENT(IN) :: MeanPeriodHours ! Averaging period in hours 
    7678      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    7779 
     
    280282               inpfiles(jj)%iobsj = -1 
    281283            ENDIF 
     284 
     285            !If SLA observations are representing a time mean then set the time 
     286            !of the obs to the end of that meaning period relative to the start of the run 
     287            IF ( ln_time_mean_sla_bkg .AND. ( TRIM( clvarsin(1) ) == 'SLA' ) ) THEN 
     288               DO ji = 1, inpfiles(jj)%nobs 
     289                  ! Only do this for obs within time window 
     290                  IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 
     291                     & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 
     292                     inpfiles(jj)%ptim(ji) = & 
     293                           & djulini(jj) + (MeanPeriodHours/24.) 
     294                  ENDIF       
     295               END DO 
     296            ENDIF    
     297 
    282298            inowin = 0 
    283299            DO ji = 1, inpfiles(jj)%nobs 
Note: See TracChangeset for help on using the changeset viewer.