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 7858 – NEMO

Changeset 7858


Ignore:
Timestamp:
2017-03-31T12:05:33+02:00 (7 years ago)
Author:
mattmartin
Message:

Commit updates to include the averaging observation operator.

Location:
branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r7837 r7858  
    4848   LOGICAL :: ln_diaobs           !: Logical switch for the obs operator 
    4949   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs 
     50   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres 
     51   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres 
     52   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres 
     53   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     54 
     55   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 
     56   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 
     57   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 
     58   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 
     59   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 
     60   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 
     61   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 
     62   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 
    5063 
    5164   INTEGER :: nn_1dint       !: Vertical interpolation method 
    52    INTEGER :: nn_2dint       !: Horizontal interpolation method 
     65   INTEGER :: nn_2dint       !: Default horizontal interpolation method 
     66   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method  
     67   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method  
     68   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method  
     69   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method  
     70  
    5371   INTEGER, DIMENSION(imaxavtypes) :: & 
    5472      & nn_profdavtypes      !: Profile data types representing a daily average 
    … …  
    6179      & nextrprof, &         !: Number of profile extra variables 
    6280      & nextrsurf            !: Number of surface extra variables 
     81   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     82      & n2dintsurf           !: Interpolation option for surface variables 
     83   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     84      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     85      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables 
     86   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
     87      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     88      & llnightav            !: Logical for calculating night-time averages 
    6389 
    6490   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
    … …  
    132158         & cn_altbiasfile        ! Altimeter bias input filename 
    133159 
    134       CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
    135          & clproffiles, &        ! Profile filenames 
    136          & clsurffiles           ! Surface filenames 
    137160 
    138161      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
    … …  
    153176      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
    154177      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
     178 
     179      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
     180      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
     181 
     182      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
     183         & clproffiles, &        ! Profile filenames 
     184         & clsurffiles           ! Surface filenames 
     185 
    155186      LOGICAL :: llvar1          ! Logical for profile variable 1 
    156187      LOGICAL :: llvar2          ! Logical for profile variable 1 
    157       LOGICAL :: llnightav       ! Logical for calculating night-time averages 
    158  
    159       REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
    160       REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
     188 
    161189      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    162190         & zglam1, &             ! Model longitudes for profile variable 1 
    … …  
    169197         & zmask2                ! Model land/sea mask associated with variable 2 
    170198 
     199 
    171200      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    172201         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
    … …  
    176205         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    177206         &            ln_sstnight,                                    & 
     207         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     208         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    178209         &            cn_profbfiles, cn_slafbfiles,                   & 
    179210         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    … …  
    183214         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    184215         &            cn_gridsearchfile, rn_gridsearchres,            & 
    185          &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     216         &            rn_dobsini, rn_dobsend,                         & 
     217         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     218         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     219         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     220         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     221         &            nn_1dint, nn_2dint,                             & 
     222         &            nn_2dint_sla, nn_2dint_sst,                     & 
     223         &            nn_2dint_sss, nn_2dint_sic,                     & 
    186224         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    187225         &            nn_profdavtypes 
    … …  
    236274      ENDIF 
    237275 
    238  
    239       !----------------------------------------------------------------------- 
    240       ! Set up list of observation types to be used 
    241       ! and the files associated with each type 
    242       !----------------------------------------------------------------------- 
    243  
    244       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    245       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 
    246          &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) 
    247  
    248       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    249          IF(lwp) WRITE(numout,cform_war) 
    250          IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
    251             &                    ' are set to .FALSE. so turning off calls to dia_obs' 
    252          nwarn = nwarn + 1 
    253          lk_diaobs = .FALSE. 
    254          RETURN 
    255       ENDIF 
    256  
    257       IF ( nproftypes > 0 ) THEN 
    258  
    259          ALLOCATE( cobstypesprof(nproftypes) ) 
    260          ALLOCATE( ifilesprof(nproftypes) ) 
    261          ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
    262  
    263          jtype = 0 
    264          IF (ln_t3d .OR. ln_s3d) THEN 
    265             jtype = jtype + 1 
    266             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
    267                &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    268          ENDIF 
    269          IF (ln_vel3d) THEN 
    270             jtype = jtype + 1 
    271             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
    272                &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    273          ENDIF 
    274  
    275       ENDIF 
    276  
    277       IF ( nsurftypes > 0 ) THEN 
    278  
    279          ALLOCATE( cobstypessurf(nsurftypes) ) 
    280          ALLOCATE( ifilessurf(nsurftypes) ) 
    281          ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
    282  
    283          jtype = 0 
    284          IF (ln_sla) THEN 
    285             jtype = jtype + 1 
    286             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
    287                &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    288          ENDIF 
    289          IF (ln_sst) THEN 
    290             jtype = jtype + 1 
    291             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
    292                &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    293          ENDIF 
    294 #if defined key_lim2 || defined key_lim3 || defined key_cice 
    295          IF (ln_sic) THEN 
    296             jtype = jtype + 1 
    297             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
    298                &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    299          ENDIF 
    300 #endif 
    301          IF (ln_sss) THEN 
    302             jtype = jtype + 1 
    303             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
    304                &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    305          ENDIF 
    306  
    307          IF (ln_logchl) THEN 
    308             jtype = jtype + 1 
    309             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', & 
    310                &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    311          ENDIF 
    312  
    313          IF (ln_spm) THEN 
    314             jtype = jtype + 1 
    315             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', & 
    316                &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    317          ENDIF 
    318  
    319          IF (ln_fco2) THEN 
    320             jtype = jtype + 1 
    321             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', & 
    322                &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    323          ENDIF 
    324  
    325          IF (ln_pco2) THEN 
    326             jtype = jtype + 1 
    327             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', & 
    328                &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    329          ENDIF 
    330  
    331       ENDIF 
    332  
    333       !Write namelist settings to stdout 
    334276      IF(lwp) THEN 
    335277         WRITE(numout,*) 
    … …  
    348290         WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2 
    349291         WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2 
    350          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    351          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
     292         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     293         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    352294         IF (ln_grid_search_lookup) & 
    353295            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
    … …  
    366308         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    367309         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    368          WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    369  
    370          IF ( nproftypes > 0 ) THEN 
    371             DO jtype = 1, nproftypes 
    372                DO jfile = 1, ifilesprof(jtype) 
    373                   WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
    374                      TRIM(clproffiles(jtype,jfile)) 
    375                END DO 
    376             END DO 
     310      ENDIF 
     311      !----------------------------------------------------------------------- 
     312      ! Set up list of observation types to be used 
     313      ! and the files associated with each type 
     314      !----------------------------------------------------------------------- 
     315 
     316      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     317      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 
     318         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) 
     319 
     320      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     321         IF(lwp) WRITE(numout,cform_war) 
     322         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     323            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     324         nwarn = nwarn + 1 
     325         lk_diaobs = .FALSE. 
     326         RETURN 
     327      ENDIF 
     328 
     329      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     330      IF ( nproftypes > 0 ) THEN 
     331 
     332         ALLOCATE( cobstypesprof(nproftypes) ) 
     333         ALLOCATE( ifilesprof(nproftypes) ) 
     334         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     335 
     336         jtype = 0 
     337         IF (ln_t3d .OR. ln_s3d) THEN 
     338            jtype = jtype + 1 
     339            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     340               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    377341         ENDIF 
    378  
    379          WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
    380          IF ( nsurftypes > 0 ) THEN 
    381             DO jtype = 1, nsurftypes 
    382                DO jfile = 1, ifilessurf(jtype) 
    383                   WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
    384                      TRIM(clsurffiles(jtype,jfile)) 
    385                END DO 
    386             END DO 
     342         IF (ln_vel3d) THEN 
     343            jtype = jtype + 1 
     344            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     345               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    387346         ENDIF 
    388          WRITE(numout,*) '~~~~~~~~~~~~' 
    389  
    390       ENDIF 
     347 
     348      ENDIF 
     349 
     350      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
     351      IF ( nsurftypes > 0 ) THEN 
     352 
     353         ALLOCATE( cobstypessurf(nsurftypes) ) 
     354         ALLOCATE( ifilessurf(nsurftypes) ) 
     355         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     356         ALLOCATE(n2dintsurf(nsurftypes)) 
     357         ALLOCATE(ravglamscl(nsurftypes)) 
     358         ALLOCATE(ravgphiscl(nsurftypes)) 
     359         ALLOCATE(lfpindegs(nsurftypes)) 
     360         ALLOCATE(llnightav(nsurftypes)) 
     361 
     362         jtype = 0 
     363         IF (ln_sla) THEN 
     364            jtype = jtype + 1 
     365            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
     366               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     367            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
     368               &                  nn_2dint, nn_2dint_sla,             & 
     369               &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
     370               &                  ln_sla_fp_indegs, .FALSE.,          & 
     371               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     372               &                  lfpindegs, llnightav ) 
     373         ENDIF 
     374         IF (ln_sst) THEN 
     375            jtype = jtype + 1 
     376            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
     377               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     378            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
     379               &                  nn_2dint, nn_2dint_sst,             & 
     380               &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
     381               &                  ln_sst_fp_indegs, ln_sstnight,      & 
     382               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     383               &                  lfpindegs, llnightav ) 
     384         ENDIF 
     385#if defined key_lim2 || defined key_lim3 || defined key_cice 
     386         IF (ln_sic) THEN 
     387            jtype = jtype + 1 
     388            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
     389               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     390            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
     391               &                  nn_2dint, nn_2dint_sic,             & 
     392               &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
     393               &                  ln_sic_fp_indegs, .FALSE.,          & 
     394               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     395               &                  lfpindegs, llnightav ) 
     396         ENDIF 
     397#endif 
     398         IF (ln_sss) THEN 
     399            jtype = jtype + 1 
     400            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
     401               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     402            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
     403               &                  nn_2dint, nn_2dint_sss,             & 
     404               &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
     405               &                  ln_sss_fp_indegs, .FALSE.,          & 
     406               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     407               &                  lfpindegs, llnightav ) 
     408         ENDIF 
     409 
     410         IF (ln_logchl) THEN 
     411            jtype = jtype + 1 
     412            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', & 
     413               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     414            CALL obs_setinterpopts( nsurftypes, jtype, 'logchl',         & 
     415               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     416               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     417               &                    lfpindegs, llnightav ) 
     418         ENDIF 
     419 
     420         IF (ln_spm) THEN 
     421            jtype = jtype + 1 
     422            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', & 
     423               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     424            CALL obs_setinterpopts( nsurftypes, jtype, 'spm   ',         & 
     425               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     426               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     427               &                    lfpindegs, llnightav ) 
     428         ENDIF 
     429 
     430         IF (ln_fco2) THEN 
     431            jtype = jtype + 1 
     432            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', & 
     433               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     434            CALL obs_setinterpopts( nsurftypes, jtype, 'fco2  ',         & 
     435               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     436               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     437               &                    lfpindegs, llnightav ) 
     438         ENDIF 
     439 
     440         IF (ln_pco2) THEN 
     441            jtype = jtype + 1 
     442            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', & 
     443               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     444            CALL obs_setinterpopts( nsurftypes, jtype, 'pco2  ',         & 
     445               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     446               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     447               &                    lfpindegs, llnightav ) 
     448         ENDIF 
     449 
     450      ENDIF 
     451 
     452      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     453 
    391454 
    392455      !----------------------------------------------------------------------- 
    … …  
    404467      ENDIF 
    405468 
    406       IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
     469      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 
    407470         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    408471            &                    ' is not available') 
    … …  
    488551            nvarssurf(jtype) = 1 
    489552            nextrsurf(jtype) = 0 
    490             llnightav = .FALSE. 
    491553            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
    492             IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 
    493554 
    494555            !Read in surface obs types 
    … …  
    496557               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    497558               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    498                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
     559               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
    499560 
    500561            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    501562 
    502563            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    503                CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
    504                IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     564               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     565               IF ( ln_altbias ) & 
     566                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    505567            ENDIF 
    506568 
    … …  
    515577               ENDIF 
    516578 
    517                CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, &  
     579               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
    518580                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
    519581 
    … …  
    623685         & zgphi1,    &            ! Model latitudes for prof variable 1 
    624686         & zgphi2                  ! Model latitudes for prof variable 2 
    625       LOGICAL :: llnightav        ! Logical for calculating night-time average 
    626687 
    627688 
    … …  
    697758            !Defaults which might be changed 
    698759            zsurfmask(:,:) = tmask(:,:,1) 
    699             llnightav = .FALSE. 
    700760 
    701761            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    702762            CASE('sst') 
    703763               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    704                llnightav = ln_sstnight 
    705764            CASE('sla') 
    706765               zsurfvar(:,:) = sshn(:,:) 
    … …  
    741800                  &           ' but no biogeochemical model appears to have been defined' ) 
    742801#endif 
    743                llnightav = .FALSE. 
    744802               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
    745803               ! Take the log10 where we can, otherwise exclude 
    … …  
    846904            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    847905               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
    848                &               nn_2dint, llnightav ) 
     906               &               n2dintsurf(jtype), llnightav(jtype),     & 
     907               &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     908               &               lfpindegs(jtype) ) 
    849909 
    850910         END DO 
    … …  
    9781038 
    9791039      IF ( nsurftypes > 0 ) & 
    980          &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     1040         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     1041         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 
    9811042 
    9821043   END SUBROUTINE dia_obs_dealloc 
    … …  
    11771238    ENDIF 
    11781239 
     1240    IF(lwp) THEN     
     1241       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1242       DO jfile = 1, ifiles(jtype) 
     1243          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1244       END DO 
     1245    ENDIF 
     1246 
    11791247    END SUBROUTINE obs_settypefiles 
    11801248 
     1249    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1250               &                  n2dint_default, n2dint_type,        & 
     1251               &                  ravglamscl_type, ravgphiscl_type,   & 
     1252               &                  lfp_indegs_type, lavnight_type,     & 
     1253               &                  n2dint, ravglamscl, ravgphiscl,     & 
     1254               &                  lfpindegs, lavnight ) 
     1255 
     1256    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1257    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1258    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1259    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1260    REAL(wp), INTENT(IN) :: & 
     1261       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1262       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1263    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1264    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1265    CHARACTER(len=6), INTENT(IN) :: ctypein  
     1266 
     1267    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1268       &                    n2dint  
     1269    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1270       &                    ravglamscl, ravgphiscl 
     1271    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1272       &                    lfpindegs, lavnight 
     1273 
     1274    lavnight(jtype) = lavnight_type 
     1275 
     1276    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 
     1277       n2dint(jtype) = n2dint_type 
     1278    ELSE 
     1279       n2dint(jtype) = n2dint_default 
     1280    ENDIF 
     1281 
     1282    ! For averaging observation footprints set options for size of footprint  
     1283    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1284       IF ( ravglamscl_type > 0._wp ) THEN 
     1285          ravglamscl(jtype) = ravglamscl_type 
     1286       ELSE 
     1287          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1288                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1289       ENDIF 
     1290 
     1291       IF ( ravgphiscl_type > 0._wp ) THEN 
     1292          ravgphiscl(jtype) = ravgphiscl_type 
     1293       ELSE 
     1294          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1295                         'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1296       ENDIF 
     1297 
     1298       lfpindegs(jtype) = lfp_indegs_type  
     1299 
     1300    ENDIF 
     1301 
     1302    ! Write out info  
     1303    IF(lwp) THEN 
     1304       IF ( n2dint(jtype) <= 4 ) THEN 
     1305          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1306             &            ' model counterparts will be interpolated horizontally' 
     1307       ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1308          WRITE(numout,*) TRIM(ctypein)// & 
     1309             &            ' model counterparts will be averaged horizontally' 
     1310          WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1311          WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1312          IF ( lfpindegs(jtype) ) THEN 
     1313              WRITE(numout,*) '             '//'    (in degrees)' 
     1314          ELSE 
     1315              WRITE(numout,*) '             '//'    (in metres)' 
     1316          ENDIF 
     1317       ENDIF 
     1318    ENDIF 
     1319 
     1320    END SUBROUTINE obs_setinterpopts 
     1321 
    11811322END MODULE diaobs 

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r7837 r7858  
    1919      & obs_int_h2d, & 
    2020      & obs_int_h2d_init 
     21   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint 
     22      & obs_avg_h2d, & 
     23      & obs_avg_h2d_init, & 
     24      & obs_max_fpsize 
    2125   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    2226      & obs_int_z1d,    & 
    … …  
    2428   USE obs_const,  ONLY :    &    ! Obs fill value 
    2529      & obfillflt 
    26    USE dom_oce,       ONLY : &    ! Model lats/lons 
    27       & glamt, & 
    28       & gphit 
     30   USE dom_oce,       ONLY : & 
     31      & glamt, glamf, & 
     32      & gphit, gphif 
    2933   USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    3034      & ctl_warn, ctl_stop 
    … …  
    667671   END SUBROUTINE obs_prof_opt 
    668672 
    669    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         & 
    670       &                    kit000, kdaystp, psurf, psurfmask, & 
    671       &                    k2dint, ldnightav ) 
     673   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     674      &                     kit000, kdaystp, psurf, psurfmask,   & 
     675      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     676      &                     lindegrees ) 
    672677 
    673678      !!----------------------------------------------------------------------- 
    … …  
    691696      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    692697      !! 
     698      !!    Two horizontal averaging schemes are also available: 
     699      !!        - weighted radial footprint        (k2dint = 5) 
     700      !!        - weighted rectangular footprint   (k2dint = 6) 
     701      !! 
    693702      !! 
    694703      !! ** Action  : 
    … …  
    697706      !!      ! 07-03 (A. Weaver) 
    698707      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     708      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    699709      !!----------------------------------------------------------------------- 
    700710 
    … …  
    718728         & psurfmask                   ! Land-sea mask 
    719729      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     730      REAL(KIND=wp), INTENT(IN) :: & 
     731         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     732         & pphiscl                     ! This is the full width (rather than half-width) 
     733      LOGICAL, INTENT(IN) :: & 
     734         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
    720735 
    721736      !! * Local declarations 
    … …  
    726741      INTEGER :: isurf 
    727742      INTEGER :: iobs 
     743      INTEGER :: imaxifp, imaxjfp 
     744      INTEGER :: imodi, imodj 
    728745      INTEGER :: idayend 
    729746      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    730          & igrdi, & 
    731          & igrdj 
     747         & igrdi,   & 
     748         & igrdj,   & 
     749         & igrdip1, & 
     750         & igrdjp1 
    732751      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    733752         & icount_night,      & 
    … …  
    737756      REAL(wp), DIMENSION(1) :: zext, zobsmask 
    738757      REAL(wp) :: zdaystp 
    739       REAL(wp), DIMENSION(2,2,1) :: & 
    740          & zweig 
    741758      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     759         & zweig,  & 
    742760         & zmask,  & 
    743761         & zsurf,  & 
    744762         & zsurfm, & 
     763         & zsurftmp, & 
    745764         & zglam,  & 
    746          & zgphi 
     765         & zgphi,  & 
     766         & zglamf, & 
     767         & zgphif 
     768 
    747769      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    748770         & zintmp,  & 
    … …  
    756778      inrc = kt - kit000 + 2 
    757779      isurf = surfdataqc%nsstp(inrc) 
     780 
     781      ! Work out the maximum footprint size for the  
     782      ! interpolation/averaging in model grid-points - has to be even. 
     783 
     784      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     785 
    758786 
    759787      IF ( ldnightav ) THEN 
    … …  
    822850 
    823851      ALLOCATE( & 
    824          & igrdi(2,2,isurf), & 
    825          & igrdj(2,2,isurf), & 
    826          & zglam(2,2,isurf), & 
    827          & zgphi(2,2,isurf), & 
    828          & zmask(2,2,isurf), & 
    829          & zsurf(2,2,isurf)  & 
     852         & zweig(imaxifp,imaxjfp,1),      & 
     853         & igrdi(imaxifp,imaxjfp,isurf), & 
     854         & igrdj(imaxifp,imaxjfp,isurf), & 
     855         & zglam(imaxifp,imaxjfp,isurf), & 
     856         & zgphi(imaxifp,imaxjfp,isurf), & 
     857         & zmask(imaxifp,imaxjfp,isurf), & 
     858         & zsurf(imaxifp,imaxjfp,isurf), & 
     859         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     860         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     861         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     862         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     863         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    830864         & ) 
    831865 
    832866      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    833867         iobs = jobs - surfdataqc%nsurfup 
    834          igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
    835          igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
    836          igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
    837          igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
    838          igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
    839          igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
    840          igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
    841          igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
     868         DO ji = 0, imaxifp 
     869            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     870             
     871            !Deal with wrap around in longitude 
     872            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     873            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     874             
     875            DO jj = 0, imaxjfp 
     876               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     877               !If model values are out of the domain to the north/south then 
     878               !set them to be the edge of the domain 
     879               IF ( imodj < 1      ) imodj = 1 
     880               IF ( imodj > jpjglo ) imodj = jpjglo 
     881 
     882               igrdip1(ji+1,jj+1,iobs) = imodi 
     883               igrdjp1(ji+1,jj+1,iobs) = imodj 
     884                
     885               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     886                  igrdi(ji,jj,iobs) = imodi 
     887                  igrdj(ji,jj,iobs) = imodj 
     888               ENDIF 
     889                
     890            END DO 
     891         END DO 
    842892      END DO 
    843893 
    844       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     894      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    845895         &                  igrdi, igrdj, glamt, zglam ) 
    846       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     896      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    847897         &                  igrdi, igrdj, gphit, zgphi ) 
    848       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     898      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    849899         &                  igrdi, igrdj, psurfmask, zmask ) 
    850       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     900      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    851901         &                  igrdi, igrdj, psurf, zsurf ) 
     902      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     903         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     904      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     905         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    852906 
    853907      ! At the end of the day get interpolated means 
    … …  
    855909 
    856910         ALLOCATE( & 
    857             & zsurfm(2,2,isurf)  & 
     911            & zsurfm(imaxifp,imaxjfp,isurf)  & 
    858912            & ) 
    859913 
    860          CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 
     914         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
    861915            &               surfdataqc%vdmean(:,:), zsurfm ) 
    862916 
    … …  
    889943         zphi = surfdataqc%rphi(jobs) 
    890944 
    891          ! Get weights to interpolate the model value to the observation point 
    892          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    893             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    894             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    895  
    896          ! Interpolate the model field to the observation point 
    897945         IF ( ldnightav .AND. idayend == 0 ) THEN 
    898946            ! Night-time averaged data 
    899             CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
     947            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    900948         ELSE 
    901             CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     949            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     950         ENDIF 
     951 
     952         IF ( k2dint <= 4 ) THEN 
     953 
     954            ! Get weights to interpolate the model value to the observation point 
     955            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     956               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     957               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     958 
     959            ! Interpolate the model value to the observation point  
     960            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     961 
     962         ELSE 
     963 
     964            ! Get weights to average the model SLA to the observation footprint 
     965            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     966               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     967               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     968               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     969               &                   lindegrees, zweig, zobsmask ) 
     970 
     971            ! Average the model SST to the observation footprint 
     972            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     973               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     974 
    902975         ENDIF 
    903976 
    … …  
    914987      ! Deallocate the data for interpolation 
    915988      DEALLOCATE( & 
     989         & zweig, & 
    916990         & igrdi, & 
    917991         & igrdj, & 
    … …  
    919993         & zgphi, & 
    920994         & zmask, & 
    921          & zsurf  & 
     995         & zsurf, & 
     996         & zsurftmp, & 
     997         & zglamf, & 
     998         & zgphif, & 
     999         & igrdip1,& 
     1000         & igrdjp1 & 
    9221001         & ) 
    9231002 
Note: See TracChangeset for help on using the changeset viewer.