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 7773 for branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Ignore:
Timestamp:
2017-03-09T13:52:43+01:00 (7 years ago)
Author:
mattmartin
Message:

Committing updates after doing the following:

  • merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
  • updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
    • generalised vertical coordinate for profile obs. This was done so that is now the default option.
    • sst bias correction implemented with the new simplified obs code.
    • included the biogeochemical obs types int he new simplified obs code.
    • included the changes to exclude obs in the boundary for limited area models
    • included other changes for the efficiency of the obs operator to remove global arrays.
Location:
branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
38 edited
36 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/ddatetoymdhms.h90

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted
    r5704 r7773  
    2727   USE obs_grid                 ! Grid searching 
    2828   USE obs_read_altbias         ! Bias treatment for altimeter 
     29   USE obs_sstbias              ! Bias correction routine for SST 
    2930   USE obs_profiles_def         ! Profile data definitions 
    3031   USE obs_surf_def             ! Surface data definitions 
    … …  
    7677   !!---------------------------------------------------------------------- 
    7778 
     79   !! * Substitutions  
     80#  include "domzgr_substitute.h90" 
    7881CONTAINS 
    7982 
    … …  
    9396      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    9497      !!        !  07-03  (K. Mogensen) General handling of profiles 
     98      !!        !  14-08  (J.While) Incorporated SST bias correction 
    9599      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    96100      !!---------------------------------------------------------------------- 
    … …  
    108112      INTEGER :: jvar            ! Counter for variables 
    109113      INTEGER :: jfile           ! Counter for files 
     114      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
    110115 
    111116      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
    112          & cn_profbfiles, &      ! T/S profile input filenames 
    113          & cn_sstfbfiles, &      ! Sea surface temperature input filenames 
    114          & cn_slafbfiles, &      ! Sea level anomaly input filenames 
    115          & cn_sicfbfiles, &      ! Seaice concentration input filenames 
    116          & cn_velfbfiles         ! Velocity profile input filenames 
     117         & cn_profbfiles,    &   ! T/S profile input filenames 
     118         & cn_sstfbfiles,    &   ! Sea surface temperature input filenames 
     119         & cn_slafbfiles,    &   ! Sea level anomaly input filenames 
     120         & cn_sicfbfiles,    &   ! Seaice concentration input filenames 
     121         & cn_velfbfiles     &   ! Velocity profile input filenames 
     122         & cn_sssfbfiles,    &   ! Sea surface salinity input filenames 
     123         & cn_logchlfbfiles, &   ! Log(Chl) input filenames 
     124         & cn_spmfbfiles,    &   ! Sediment input filenames 
     125         & cn_fco2fbfiles,   &   ! fco2 input filenames 
     126         & cn_pco2fbfiles,   &   ! pco2 input filenames 
     127         & cn_sstbiasfiles       ! SST bias input filenames 
     128 
    117129      CHARACTER(LEN=128) :: & 
    118130         & cn_altbiasfile        ! Altimeter bias input filename 
     131 
    119132      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
    120133         & clproffiles, &        ! Profile filenames 
    … …  
    126139      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
    127140      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     141      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
    128142      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     143      LOGICAL :: ln_logchl       ! Logical switch for log(Chl) obs 
     144      LOGICAL :: ln_spm          ! Logical switch for sediment obs 
     145      LOGICAL :: ln_fco2         ! Logical switch for fco2 obs 
     146      LOGICAL :: ln_pco2         ! Logical switch for pco2 obs 
    129147      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
    130148      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
     149      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
    131150      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
    132151      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     152      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
    133153      LOGICAL :: llvar1          ! Logical for profile variable 1 
    134154      LOGICAL :: llvar2          ! Logical for profile variable 1 
    … …  
    148168 
    149169      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    150          &            ln_sst, ln_sic, ln_vel3d,                       & 
    151          &            ln_altbias, ln_nea, ln_grid_global,             & 
    152          &            ln_grid_search_lookup,                          & 
    153          &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     170         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     171         &            ln_logchl, ln_spm, ln_fco2, ln_pco2,            & 
     172         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     173         &            ln_grid_global, ln_grid_search_lookup,          & 
     174         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     175 
     176ln_sstnight,              & 
    154177         &            cn_profbfiles, cn_slafbfiles,                   & 
    155178         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    156          &            cn_velfbfiles, cn_altbiasfile,                  & 
     179         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     180         &            cn_logchlfbfiles, cn_spmfbfiles,                & 
     181         &            cn_fco2fbfiles, cn_pco2fbfiles,                 & 
     182         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    157183         &            cn_gridsearchfile, rn_gridsearchres,            & 
    158184         &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
    … …  
    172198 
    173199      ! Some namelist arrays need initialising 
    174       cn_profbfiles(:) = '' 
    175       cn_slafbfiles(:) = '' 
    176       cn_sstfbfiles(:) = '' 
    177       cn_sicfbfiles(:) = '' 
    178       cn_velfbfiles(:) = '' 
    179       nn_profdavtypes(:) = -1 
     200      cn_profbfiles(:)    = '' 
     201      cn_slafbfiles(:)    = '' 
     202      cn_sstfbfiles(:)    = '' 
     203      cn_sicfbfiles(:)    = '' 
     204      cn_velfbfiles(:)    = '' 
     205      cn_sssfbfiles(:)    = '' 
     206      cn_logchlfbfiles(:) = '' 
     207      cn_spmfbfiles(:)    = '' 
     208      cn_fco2fbfiles(:)   = '' 
     209      cn_pco2fbfiles(:)   = '' 
     210      cn_sstbiasfiles(:)  = '' 
     211      nn_profdavtypes(:)  = -1 
    180212 
    181213      CALL ini_date( rn_dobsini ) 
    … …  
    204236 
    205237      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    206       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
     238      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 
     239         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) 
    207240 
    208241      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    … …  
    285318         ENDIF 
    286319#endif 
     320         IF (ln_sss) THEN 
     321            jtype = jtype + 1 
     322            clsurffiles(jtype,:) = cn_sssfbfiles(:) 
     323            cobstypessurf(jtype) = 'sss   ' 
     324            ifilessurf(jtype) = 0 
     325            DO jfile = 1, jpmaxnfiles 
     326               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     327                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     328            END DO 
     329         ENDIF 
     330 
     331         IF (ln_logchl) THEN 
     332            jtype = jtype + 1 
     333            clsurffiles(jtype,:) = cn_logchlfbfiles(:) 
     334            cobstypessurf(jtype) = 'logchl' 
     335            ifilessurf(jtype) = 0 
     336            DO jfile = 1, jpmaxnfiles 
     337               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     338                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     339            END DO 
     340         ENDIF 
     341 
     342         IF (ln_spm) THEN 
     343            jtype = jtype + 1 
     344            clsurffiles(jtype,:) = cn_spmfbfiles(:) 
     345            cobstypessurf(jtype) = 'spm   ' 
     346            ifilessurf(jtype) = 0 
     347            DO jfile = 1, jpmaxnfiles 
     348               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     349                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     350            END DO 
     351         ENDIF 
     352 
     353         IF (ln_fco2) THEN 
     354            jtype = jtype + 1 
     355            clsurffiles(jtype,:) = cn_fco2fbfiles(:) 
     356            cobstypessurf(jtype) = 'fco2  ' 
     357            ifilessurf(jtype) = 0 
     358            DO jfile = 1, jpmaxnfiles 
     359               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     360                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     361            END DO 
     362         ENDIF 
     363 
     364         IF (ln_pco2) THEN 
     365            jtype = jtype + 1 
     366            clsurffiles(jtype,:) = cn_pco2fbfiles(:) 
     367            cobstypessurf(jtype) = 'pco2  ' 
     368            ifilessurf(jtype) = 0 
     369            DO jfile = 1, jpmaxnfiles 
     370               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
     371                  ifilessurf(jtype) = ifilessurf(jtype) + 1 
     372            END DO 
     373         ENDIF 
    287374 
    288375      ENDIF 
    … …  
    300387         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    301388         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
     389         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     390         WRITE(numout,*) '             Logical switch for log(Chl) observations              ln_logchl = ', ln_logchl 
     391         WRITE(numout,*) '             Logical switch for SPM observations                      ln_spm = ', ln_spm 
     392         WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2 
     393         WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2 
    302394         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    303395         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
    … …  
    309401         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
    310402         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     403         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    311404         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    312405         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    313406         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    314407         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     408         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    315409         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    316410         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    … …  
    418512               &               jpi, jpj, jpk, & 
    419513               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
    420                &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     514               &               ln_nea, ln_bound_reject, & 
     515               &               kdailyavtypes = nn_profdavtypes ) 
    421516 
    422517         END DO 
    … …  
    447542               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    448543 
    449             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     544            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    450545 
    451546            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    … …  
    453548               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
    454549            ENDIF 
     550 
     551            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
     552               jnumsstbias = 0 
     553               DO jfile = 1, jpmaxnfiles 
     554                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     555                  jnumsstbias = jnumsstbias + 1 
     556               END DO 
     557               IF ( jnumsstbias == 0 ) THEN 
     558                  CALL ctl_stop("ln_sstbias set,"// &  
     559                     &          "  but no bias files to read in")     
     560               ENDIF 
     561 
     562               CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, &  
     563                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
    455564 
    456565         END DO 
    … …  
    507616         & frld 
    508617#endif 
     618#if defined key_hadocc 
     619      USE trc, ONLY :  &                ! HadOCC chlorophyll, fCO2 and pCO2 
     620         & HADOCC_CHL, & 
     621         & HADOCC_FCO2, & 
     622         & HADOCC_PCO2, & 
     623         & HADOCC_FILL_FLT 
     624#elif defined key_medusa && defined key_foam_medusa 
     625      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2 
     626         & MEDUSA_CHL, & 
     627         & MEDUSA_FCO2, & 
     628         & MEDUSA_PCO2, & 
     629         & MEDUSA_FILL_FLT 
     630#elif defined key_fabm 
     631      USE fabm 
     632      USE par_fabm 
     633#endif 
     634#if defined key_spm 
     635      USE par_spm, ONLY: &              ! ERSEM/SPM sediments 
     636         & jp_spm 
     637      USE trc, ONLY :  & 
     638         & trn 
     639#endif 
     640 
    509641      IMPLICIT NONE 
    510642 
    … …  
    523655         & zprofmask2              ! Mask associated with zprofvar2 
    524656      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    525          & zsurfvar                ! Model values equivalent to surface ob. 
     657         & zsurfvar, &             ! Model values equivalent to surface ob. 
     658         & zsurfmask               ! Mask associated with surface variable 
    526659      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    527660         & zglam1,    &            ! Model longitudes for prof variable 1 
    … …  
    540673      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
    541674      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     675      CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    542676      CALL wrk_alloc( jpi, jpj, zglam1 ) 
    543677      CALL wrk_alloc( jpi, jpj, zglam2 ) 
    … …  
    608742         DO jtype = 1, nsurftypes 
    609743 
     744            !Defaults which might be changed 
     745            zsurfmask(:,:) = tmask(:,:,1) 
     746            llnightav = .FALSE. 
     747 
    610748            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    611749            CASE('sst') 
    … …  
    614752            CASE('sla') 
    615753               zsurfvar(:,:) = sshn(:,:) 
    616                llnightav = .FALSE. 
     754            CASE('sss') 
     755               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    617756#if defined key_lim2 || defined key_lim3 
    618757            CASE('sic') 
    … …  
    630769                  zsurfvar(:,:) = 1._wp - frld(:,:) 
    631770               ENDIF 
    632  
     771#endif 
     772            CASE('logchl') 
     773#if defined key_hadocc 
     774               zsurfvar(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     775#elif defined key_medusa && defined key_foam_medusa 
     776               zsurfvar(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     777#elif defined key_fabm 
     778               chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot) 
     779               zsurfvar(:,:) = chl_3d(:,:,1) 
     780#else 
     781               CALL ctl_stop( ' Trying to run logchl observation operator', & 
     782                  &           ' but no biogeochemical model appears to have been defined' ) 
     783#endif 
    633784               llnightav = .FALSE. 
    634 #endif 
     785               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
     786               ! Take the log10 where we can, otherwise exclude 
     787               tiny = 1.0e-20 
     788               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
     789                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
     790               ELSEWHERE 
     791                  zsurfvar(:,:)  = obfillflt 
     792                  zsurfmask(:,:) = 0 
     793               END WHERE 
     794            CASE('spm') 
     795#if defined key_spm 
     796               zsurfvar(:,:) = 0.0 
     797               DO jn = 1, jp_spm 
     798                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     799               END DO 
     800#else 
     801               CALL ctl_stop( ' Trying to run spm observation operator', & 
     802                  &           ' but no spm model appears to have been defined' ) 
     803#endif 
     804            CASE('fco2') 
     805#if defined key_hadocc 
     806               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     807               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     808                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     809                  zsurfvar(:,:) = obfillflt 
     810                  zsurfmask(:,:) = 0 
     811                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
     812                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     813                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
     814               ENDIF 
     815#elif defined key_medusa && defined key_foam_medusa 
     816               zsurfmask(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA 
     817               IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) .AND. & 
     818                  & ( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     819                  zsurfvar(:,:) = obfillflt 
     820                  zsurfmask(:,:) = 0 
     821                  CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', & 
     822                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     823                     &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' ) 
     824               ENDIF 
     825#elif defined key_fabm 
     826               ! First, get pCO2 from FABM 
     827               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     828               zsurfvar(:,:) = pco2_3d(:,:,1) 
     829               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
     830               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
     831               ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
     832               ! and 
     833               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
     834               ! Marine Chemistry, 2: 203-215. 
     835               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
     836               ! not explicitly included - atmospheric pressure is not necessarily available so this is 
     837               ! the best assumption. 
     838               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
     839               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
     840               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
     841               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
     842               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
     843                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     844                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     845                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     846                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     847                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     848#else 
     849               CALL ctl_stop( ' Trying to run fco2 observation operator', & 
     850                  &           ' but no biogeochemical model appears to have been defined' ) 
     851#endif 
     852            CASE('pco2') 
     853#if defined key_hadocc 
     854               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     855               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     856                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     857                  zsurfvar(:,:) = obfillflt 
     858                  zsurfmask(:,:) = 0 
     859                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
     860                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     861                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
     862               ENDIF 
     863#elif defined key_medusa && defined key_foam_medusa 
     864               zsurfvar(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA 
     865               IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) .AND. & 
     866                  & ( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     867                  zsurfvar(:,:) = obfillflt 
     868                  zsurfmask(:,:) = 0 
     869                  CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', & 
     870                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     871                     &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' ) 
     872               ENDIF 
     873#elif defined key_fabm 
     874               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     875               zsurfvar(:,:) = pco2_3d(:,:,1) 
     876#else 
     877               CALL ctl_stop( ' Trying to run pCO2 observation operator', & 
     878                  &           ' but no biogeochemical model appears to have been defined' ) 
     879#endif 
     880 
    635881            END SELECT 
    636882 
    637883            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    638                &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
     884               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
    639885               &               nn_2dint, llnightav ) 
    640886 
    … …  
    648894      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
    649895      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     896      CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    650897      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    651898      CALL wrk_dealloc( jpi, jpj, zglam2 ) 

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/find_obs_proc.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/greg2jul.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/grt_cir_dis.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/grt_cir_dis_saa.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/jul2greg.h90

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/linquad.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/maxdist.h90

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_conv_functions.h90

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90

    r2358 r7773  
    325325         CALL obs_mpp_max_integer( kobsj, kobs ) 
    326326      ELSE 
    327          CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs ) 
     327         CALL obs_mpp_find_obs_proc( kproc,kobs ) 
    328328      ENDIF 
    329329 

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

    • Property svn:keywords deleted
    r5682 r7773  
    8787   !!---------------------------------------------------------------------- 
    8888   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    89    !! $Id$ 
     89   !! $Id: obs_grid.F90 5682 2015-08-12 15:46:45Z mattmartin $ 
    9090   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    9191   !!---------------------------------------------------------------------- 
    … …  
    613613         CALL obs_mpp_max_integer( kobsj, kobs ) 
    614614      ELSE 
    615          CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs ) 
     615         CALL obs_mpp_find_obs_proc( kproc, kobs ) 
    616616      ENDIF 
    617617 

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted
    r5682 r7773  
    77   !!             -   ! 2006-05  (K. Mogensen)  Reformatted 
    88   !!             -   ! 2008-01  (K. Mogensen)  add mpp_global_max 
     9   !!            3.6  ! 2015-01  (J. Waters) obs_mpp_find_obs_proc  
     10   !!                            rewritten to avoid global arrays 
    911   !!---------------------------------------------------------------------- 
    1012#  define mpivar mpi_double_precision 
    … …  
    1214   !! obs_mpp_bcast_integer : Broadcast an integer array from a processor to all processors 
    1315   !! obs_mpp_max_integer   : Find maximum on all processors of each value in an integer on all processors 
    14    !! obs_mpp_find_obs_proc : Find processors which should hold the observations 
     16   !! obs_mpp_find_obs_proc : Find processors which should hold the observations, avoiding global arrays 
    1517   !! obs_mpp_sum_integers  : Sum an integer array from all processors 
    1618   !! obs_mpp_sum_integer   : Sum an integer from all processors 
    … …  
    3739   !!---------------------------------------------------------------------- 
    3840   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    39    !! $Id$ 
     41   !! $Id: obs_mpp.F90 5682 2015-08-12 15:46:45Z mattmartin $ 
    4042   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4143   !!---------------------------------------------------------------------- 
    … …  
    9698      ! 
    9799      INTEGER :: ierr  
    98       INTEGER, DIMENSION(kno) ::   ivals 
    99       ! 
    100 INCLUDE 'mpif.h' 
    101       !!---------------------------------------------------------------------- 
     100      INTEGER, DIMENSION(:), ALLOCATABLE ::   ivals 
     101      ! 
     102INCLUDE 'mpif.h' 
     103      !!---------------------------------------------------------------------- 
     104 
     105      ALLOCATE( ivals(kno) ) 
    102106 
    103107      ! Call the MPI library to find the maximum across processors 
    … …  
    105109         &                mpi_max, mpi_comm_opa, ierr ) 
    106110      kvals(:) = ivals(:) 
     111 
     112      DEALLOCATE( ivals ) 
    107113#else 
    108114      ! no MPI: empty routine 
    … …  
    111117 
    112118 
    113    SUBROUTINE obs_mpp_find_obs_proc( kobsp, kobsi, kobsj, kno ) 
    114       !!---------------------------------------------------------------------- 
    115       !!               ***  ROUTINE obs_mpp_find_obs_proc *** 
    116       !!           
     119   SUBROUTINE obs_mpp_find_obs_proc( kobsp,kno ) 
     120      !!---------------------------------------------------------------------- 
     121      !!               ***  ROUTINE obs_mpp_find_obs_proc  *** 
     122      !!          
    117123      !! ** Purpose : From the array kobsp containing the results of the 
    118124      !!              grid search on each processor the processor return a 
    119125      !!              decision of which processors should hold the observation. 
    120126      !! 
    121       !! ** Method  : A temporary 2D array holding all the decisions is 
    122       !!              constructed using mpi_allgather on each processor. 
    123       !!              If more than one processor has found the observation 
    124       !!              with the observation in the inner domain gets it 
    125       !! 
    126       !! ** Action  : This does only work for MPI.  
     127      !! ** Method  : Synchronize the processor number for each obs using 
     128      !!              obs_mpp_max_integer. If an observation exists on two  
     129      !!              processors it will be allocated to the lower numbered 
     130      !!              processor. 
     131      !! 
     132      !! ** Action  : This does only work for MPI. 
    127133      !!              It does not work for SHMEM. 
    128134      !! 
    … …  
    130136      !!---------------------------------------------------------------------- 
    131137      INTEGER                , INTENT(in   ) ::   kno 
    132       INTEGER, DIMENSION(kno), INTENT(in   ) ::   kobsi, kobsj 
    133138      INTEGER, DIMENSION(kno), INTENT(inout) ::   kobsp 
    134139      ! 
    135140#if defined key_mpp_mpi 
    136141      ! 
    137       INTEGER :: ji 
    138       INTEGER :: jj 
    139       INTEGER :: size 
    140       INTEGER :: ierr 
    141       INTEGER :: iobsip 
    142       INTEGER :: iobsjp 
    143       INTEGER :: num_sus_obs 
    144       INTEGER, DIMENSION(kno) ::   iobsig, iobsjg 
    145       INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iobsp, iobsi, iobsj 
    146       !! 
    147 INCLUDE 'mpif.h' 
    148       !!---------------------------------------------------------------------- 
    149  
    150       !----------------------------------------------------------------------- 
    151       ! Call the MPI library to find the maximum accross processors 
    152       !----------------------------------------------------------------------- 
    153       CALL mpi_comm_size( mpi_comm_opa, size, ierr ) 
    154       !----------------------------------------------------------------------- 
    155       ! Convert local grids points to global grid points 
    156       !----------------------------------------------------------------------- 
     142      ! 
     143      INTEGER :: ji, isum 
     144      INTEGER, DIMENSION(:), ALLOCATABLE ::   iobsp 
     145      !! 
     146      !! 
     147 
     148      ALLOCATE( iobsp(kno) ) 
     149 
     150      iobsp(:)=kobsp(:) 
     151 
     152      WHERE( iobsp(:) == -1 ) 
     153         iobsp(:) = 9999999 
     154      END WHERE 
     155 
     156      iobsp(:)=-1*iobsp(:) 
     157 
     158      CALL obs_mpp_max_integer( iobsp, kno ) 
     159 
     160      kobsp(:)=-1*iobsp(:) 
     161 
     162      isum=0 
    157163      DO ji = 1, kno 
    158          IF ( ( kobsi(ji) >= 1 ) .AND. ( kobsi(ji) <= jpi ) .AND. & 
    159             & ( kobsj(ji) >= 1 ) .AND. ( kobsj(ji) <= jpj ) ) THEN 
    160             iobsig(ji) = mig( kobsi(ji) ) 
    161             iobsjg(ji) = mjg( kobsj(ji) ) 
    162          ELSE 
    163             iobsig(ji) = -1 
    164             iobsjg(ji) = -1 
     164         IF ( kobsp(ji) == 9999999 ) THEN 
     165            isum=isum+1 
     166            kobsp(ji)=-1 
    165167         ENDIF 
    166       END DO 
    167       !----------------------------------------------------------------------- 
    168       ! Get the decisions from all processors 
    169       !----------------------------------------------------------------------- 
    170       ALLOCATE( iobsp(kno,size) ) 
    171       ALLOCATE( iobsi(kno,size) ) 
    172       ALLOCATE( iobsj(kno,size) ) 
    173       CALL mpi_allgather( kobsp, kno, mpi_integer, & 
    174          &                iobsp, kno, mpi_integer, & 
    175          &                mpi_comm_opa, ierr ) 
    176       CALL mpi_allgather( iobsig, kno, mpi_integer, & 
    177          &                iobsi, kno, mpi_integer, & 
    178          &                mpi_comm_opa, ierr ) 
    179       CALL mpi_allgather( iobsjg, kno, mpi_integer, & 
    180          &                iobsj, kno, mpi_integer, & 
    181          &                mpi_comm_opa, ierr ) 
    182  
    183       !----------------------------------------------------------------------- 
    184       ! Find the processor with observations from the lowest processor  
    185       ! number among processors holding the observation. 
    186       !----------------------------------------------------------------------- 
    187       kobsp(:) = -1 
    188       num_sus_obs = 0 
    189       DO ji = 1, kno 
    190          DO jj = 1, size 
    191             IF ( ( kobsp(ji) == -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 
    192                kobsp(ji) = iobsp(ji,jj) 
    193                iobsip = iobsi(ji,jj) 
    194                iobsjp = iobsj(ji,jj) 
    195             ENDIF 
    196             IF ( ( kobsp(ji) /= -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 
    197                IF ( ( iobsip /= iobsi(ji,jj) ) .OR. & 
    198                   & ( iobsjp /= iobsj(ji,jj) ) ) THEN 
    199                   IF ( ( kobsp(ji) < 1000000 ) .AND. & 
    200                      & ( iobsp(ji,jj) < 1000000 ) ) THEN 
    201                      num_sus_obs=num_sus_obs+1 
    202                   ENDIF 
    203                ENDIF 
    204                IF ( mppmap(iobsip,iobsjp) /= ( kobsp(ji)+1 ) ) THEN 
    205                   IF ( ( iobsi(ji,jj) /= -1 ) .AND. & 
    206                      & ( iobsj(ji,jj) /= -1 ) ) THEN 
    207                      IF ((mppmap(iobsi(ji,jj),iobsj(ji,jj)) == (iobsp(ji,jj)+1))& 
    208                         & .OR. ( iobsp(ji,jj) < kobsp(ji) ) ) THEN 
    209                         kobsp(ji) = iobsp(ji,jj) 
    210                         iobsip = iobsi(ji,jj) 
    211                         iobsjp = iobsj(ji,jj) 
    212                      ENDIF 
    213                   ENDIF 
    214                ENDIF 
    215             ENDIF 
    216          END DO 
    217       END DO 
    218       IF (lwp) WRITE(numout,*) 'Number of suspicious observations: ',num_sus_obs 
    219  
    220       DEALLOCATE( iobsj ) 
    221       DEALLOCATE( iobsi ) 
     168      ENDDO 
     169 
     170 
     171      IF ( isum > 0 ) THEN 
     172         IF (lwp) WRITE(numout,*) isum, ' observations failed the grid search.' 
     173         IF (lwp) WRITE(numout,*)'If ln_grid_search_lookup=.TRUE., try reducing grid_search_res' 
     174      ENDIF 
     175 
    222176      DEALLOCATE( iobsp ) 
     177 
    223178#else 
    224179      ! no MPI: empty routine 
    225 #endif 
    226       ! 
     180#endif      
     181       
    227182   END SUBROUTINE obs_mpp_find_obs_proc 
    228183 

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

    • Property svn:keywords deleted
    r5704 r7773  
    4949   !!---------------------------------------------------------------------- 
    5050 
     51   !! * Substitutions  
     52#  include "domzgr_substitute.h90"  
    5153CONTAINS 
    5254 
    5355   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    5456      &                     kit000, kdaystp,                      & 
    55       &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
     57      &                     pvar1, pvar2, pgdept, pgdepw, 
     58      &                     pmask1, pmask2, & 
    5659      &                     plam1, plam2, pphi1, pphi2,           & 
    5760      &                     k1dint, k2dint, kdailyavtypes ) 
    … …  
    104107      !!      ! 07-03 (K. Mogensen) General handling of profiles 
    105108      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
     109      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    106110      !!----------------------------------------------------------------------- 
    107111 
    … …  
    133137         & pphi1,    &               ! Model latitudes for variable 1 
    134138         & pphi2                     ! Model latitudes for variable 2 
    135       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    136          & pgdept                    ! Model array of depth levels 
     139      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     140         & pgdept, &                 ! Model array of depth T levels  
     141         & pgdepw                    ! Model array of depth W levels  
    137142      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    138143         & kdailyavtypes             ! Types for daily averages 
    … …  
    164169         & zobsk,    & 
    165170         & zobs2k 
    166       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     171      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    167172         & zweig1, & 
    168          & zweig2 
     173         & zweig2, & 
     174         & zweig 
    169175      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    170176         & zmask1, & 
    171177         & zmask2, & 
    172          & zint1, & 
    173          & zint2, & 
    174          & zinm1, & 
    175          & zinm2 
     178         & zint1,  & 
     179         & zint2,  & 
     180         & zinm1,  & 
     181         & zinm2,  & 
     182         & zgdept, &  
     183         & zgdepw 
    176184      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    177185         & zglam1, & 
    … …  
    179187         & zgphi1, & 
    180188         & zgphi2 
     189      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2    
     190      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
     191 
    181192      LOGICAL :: ld_dailyav 
    182193 
    … …  
    259270         & zmask1(2,2,kpk,ipro),  & 
    260271         & zmask2(2,2,kpk,ipro),  & 
    261          & zint1(2,2,kpk,ipro),  & 
    262          & zint2(2,2,kpk,ipro)   & 
     272         & zint1(2,2,kpk,ipro),   & 
     273         & zint2(2,2,kpk,ipro),   & 
     274         & zgdept(2,2,kpk,ipro),  &  
     275         & zgdepw(2,2,kpk,ipro)   &  
    263276         & ) 
    264277 
    … …  
    283296      END DO 
    284297 
     298      ! Initialise depth arrays 
     299      zgdept(:,:,:,:) = 0.0 
     300      zgdepw(:,:,:,:) = 0.0 
     301 
    285302      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    286303      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
    … …  
    293310      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    294311 
     312      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     313      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     314 
    295315      ! At the end of the day also get interpolated means 
    296316      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    … …  
    307327 
    308328      ENDIF 
     329 
     330      ! Return if no observations to process  
     331      ! Has to be done after comm commands to ensure processors  
     332      ! stay in sync  
     333      IF ( ipro == 0 ) RETURN  
    309334 
    310335      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    … …  
    332357         zphi = prodatqc%rphi(jobs) 
    333358 
    334          ! Horizontal weights and vertical mask 
    335  
     359         ! Horizontal weights  
     360         ! Masked values are calculated later.   
    336361         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    337362 
    338             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     363            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    339364               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    340                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     365               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    341366 
    342367         ENDIF 
    … …  
    344369         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    345370 
    346             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     371            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    347372               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    348                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     373               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    349374  
    350375         ENDIF 
    … …  
    358383               IF ( idayend == 0 )  THEN 
    359384                  ! Daily averaged data 
    360                   CALL obs_int_h2d( kpk, kpk,      & 
    361                      &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    362  
    363                ENDIF 
    364  
    365             ELSE  
    366  
    367                ! Point data 
    368                CALL obs_int_h2d( kpk, kpk,      & 
    369                   &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    370  
    371             ENDIF 
    372  
    373             !------------------------------------------------------------- 
    374             ! Compute vertical second-derivative of the interpolating  
    375             ! polynomial at obs points 
    376             !------------------------------------------------------------- 
    377  
    378             IF ( k1dint == 1 ) THEN 
    379                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    380                   &                  pgdept, zobsmask1 ) 
    381             ENDIF 
    382  
    383             !----------------------------------------------------------------- 
    384             !  Vertical interpolation to the observation point 
    385             !----------------------------------------------------------------- 
    386             ista = prodatqc%npvsta(jobs,1) 
    387             iend = prodatqc%npvend(jobs,1) 
    388             CALL obs_int_z1d( kpk,                & 
    389                & prodatqc%var(1)%mvk(ista:iend),  & 
    390                & k1dint, iend - ista + 1,         & 
    391                & prodatqc%var(1)%vdep(ista:iend), & 
    392                & zobsk, zobs2k,                   & 
    393                & prodatqc%var(1)%vmod(ista:iend), & 
    394                & pgdept, zobsmask1 ) 
    395  
    396          ENDIF 
    397  
     385 
     386                  ! vertically interpolate all 4 corners  
     387                  ista = prodatqc%npvsta(jobs,1)  
     388                  iend = prodatqc%npvend(jobs,1)  
     389                  inum_obs = iend - ista + 1  
     390                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     391 
     392                  DO iin=1,2  
     393                     DO ijn=1,2  
     394 
     395                        IF ( k1dint == 1 ) THEN  
     396                           CALL obs_int_z1d_spl( kpk, &  
     397                              &     zinm1(iin,ijn,:,iobs), &  
     398                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     399                              &     zmask1(iin,ijn,:,iobs))  
     400                        ENDIF  
     401        
     402                        CALL obs_level_search(kpk, &  
     403                           &    zgdept(iin,ijn,:,iobs), &  
     404                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     405                           &    iv_indic)  
     406 
     407                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     408                           &    prodatqc%var(1)%vdep(ista:iend), &  
     409                           &    zinm1(iin,ijn,:,iobs), &  
     410                           &    zobs2k, interp_corner(iin,ijn,:), &  
     411                           &    zgdept(iin,ijn,:,iobs), &  
     412                           &    zmask1(iin,ijn,:,iobs))  
     413        
     414                     ENDDO  
     415                  ENDDO  
     416 
     417               ENDIF !idayend 
     418 
     419            ELSE    
     420 
     421               ! Point data  
     422      
     423               ! vertically interpolate all 4 corners  
     424               ista = prodatqc%npvsta(jobs,1)  
     425               iend = prodatqc%npvend(jobs,1)  
     426               inum_obs = iend - ista + 1  
     427               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     428               DO iin=1,2   
     429                  DO ijn=1,2  
     430                     
     431                     IF ( k1dint == 1 ) THEN  
     432                        CALL obs_int_z1d_spl( kpk, &  
     433                           &    zint1(iin,ijn,:,iobs),&  
     434                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     435                           &    zmask1(iin,ijn,:,iobs))  
     436   
     437                     ENDIF  
     438        
     439                     CALL obs_level_search(kpk, &  
     440                         &        zgdept(iin,ijn,:,iobs),&  
     441                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     442                         &        iv_indic)  
     443 
     444                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     445                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     446                         &          zint1(iin,ijn,:,iobs),            &  
     447                         &          zobs2k,interp_corner(iin,ijn,:), &  
     448                         &          zgdept(iin,ijn,:,iobs),         &  
     449                         &          zmask1(iin,ijn,:,iobs) )       
     450          
     451                  ENDDO  
     452               ENDDO  
     453              
     454            ENDIF  
     455 
     456            !-------------------------------------------------------------  
     457            ! Compute the horizontal interpolation for every profile level  
     458            !-------------------------------------------------------------  
     459              
     460            DO ikn=1,inum_obs  
     461               iend=ista+ikn-1 
     462                   
     463               zweig(:,:,1) = 0._wp  
     464    
     465               ! This code forces the horizontal weights to be   
     466               ! zero IF the observation is below the bottom of the   
     467               ! corners of the interpolation nodes, Or if it is in   
     468               ! the mask. This is important for observations near   
     469               ! steep bathymetry  
     470               DO iin=1,2  
     471                  DO ijn=1,2  
     472      
     473                     depth_loop1: DO ik=kpk,2,-1  
     474                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     475                             
     476                           zweig(iin,ijn,1) = &   
     477                              & zweig1(iin,ijn,1) * &  
     478                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     479                              &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     480                             
     481                           EXIT depth_loop1  
     482 
     483                        ENDIF  
     484 
     485                     ENDDO depth_loop1  
     486      
     487                  ENDDO  
     488               ENDDO  
     489    
     490               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     491                  &              prodatqc%var(1)%vmod(iend:iend) )  
     492 
     493                  ! Set QC flag for any observations found below the bottom 
     494                  ! needed as the check here is more strict than that in obs_prep 
     495               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     496  
     497            ENDDO  
     498  
     499            DEALLOCATE(interp_corner,iv_indic)  
     500           
     501         ENDIF  
     502 
     503         ! For the second variable 
    398504         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    399505 
    … …  
    403509 
    404510               IF ( idayend == 0 )  THEN 
    405  
    406511                  ! Daily averaged data 
    407                   CALL obs_int_h2d( kpk, kpk,      & 
    408                      &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    409  
    410                ENDIF 
    411  
    412             ELSE 
    413  
    414                ! Point data 
    415                CALL obs_int_h2d( kpk, kpk,      & 
    416                   &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    417  
    418             ENDIF 
    419  
    420  
    421             !------------------------------------------------------------- 
    422             ! Compute vertical second-derivative of the interpolating  
    423             ! polynomial at obs points 
    424             !------------------------------------------------------------- 
    425  
    426             IF ( k1dint == 1 ) THEN 
    427                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    428                   &                  pgdept, zobsmask2 ) 
    429             ENDIF 
    430  
    431             !---------------------------------------------------------------- 
    432             !  Vertical interpolation to the observation point 
    433             !---------------------------------------------------------------- 
    434             ista = prodatqc%npvsta(jobs,2) 
    435             iend = prodatqc%npvend(jobs,2) 
    436             CALL obs_int_z1d( kpk, & 
    437                & prodatqc%var(2)%mvk(ista:iend),& 
    438                & k1dint, iend - ista + 1, & 
    439                & prodatqc%var(2)%vdep(ista:iend),& 
    440                & zobsk, zobs2k, & 
    441                & prodatqc%var(2)%vmod(ista:iend),& 
    442                & pgdept, zobsmask2 ) 
    443  
    444          ENDIF 
    445  
    446       END DO 
     512 
     513                  ! vertically interpolate all 4 corners  
     514                  ista = prodatqc%npvsta(jobs,2)  
     515                  iend = prodatqc%npvend(jobs,2)  
     516                  inum_obs = iend - ista + 1  
     517                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     518 
     519                  DO iin=1,2  
     520                     DO ijn=1,2  
     521 
     522                        IF ( k1dint == 1 ) THEN  
     523                           CALL obs_int_z1d_spl( kpk, &  
     524                              &     zinm2(iin,ijn,:,iobs), &  
     525                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     526                              &     zmask2(iin,ijn,:,iobs))  
     527                        ENDIF  
     528        
     529                        CALL obs_level_search(kpk, &  
     530                           &    zgdept(iin,ijn,:,iobs), &  
     531                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     532                           &    iv_indic)  
     533 
     534                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     535                           &    prodatqc%var(2)%vdep(ista:iend), &  
     536                           &    zinm2(iin,ijn,:,iobs), &  
     537                           &    zobs2k, interp_corner(iin,ijn,:), &  
     538                           &    zgdept(iin,ijn,:,iobs), &  
     539                           &    zmask2(iin,ijn,:,iobs))  
     540        
     541                     ENDDO  
     542                  ENDDO  
     543 
     544               ENDIF !idayend 
     545 
     546            ELSE    
     547 
     548               ! Point data  
     549      
     550               ! vertically interpolate all 4 corners  
     551               ista = prodatqc%npvsta(jobs,2)  
     552               iend = prodatqc%npvend(jobs,2)  
     553               inum_obs = iend - ista + 1  
     554               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     555               DO iin=1,2   
     556                  DO ijn=1,2  
     557                     
     558                     IF ( k1dint == 1 ) THEN  
     559                        CALL obs_int_z1d_spl( kpk, &  
     560                           &    zint2(iin,ijn,:,iobs),&  
     561                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     562                           &    zmask2(iin,ijn,:,iobs))  
     563   
     564                     ENDIF  
     565        
     566                     CALL obs_level_search(kpk, &  
     567                         &        zgdept(iin,ijn,:,iobs),&  
     568                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     569                         &        iv_indic)  
     570 
     571                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     572                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     573                         &          zint2(iin,ijn,:,iobs),            &  
     574                         &          zobs2k,interp_corner(iin,ijn,:), &  
     575                         &          zgdept(iin,ijn,:,iobs),         &  
     576                         &          zmask2(iin,ijn,:,iobs) )       
     577          
     578                  ENDDO  
     579               ENDDO  
     580              
     581            ENDIF  
     582 
     583            !-------------------------------------------------------------  
     584            ! Compute the horizontal interpolation for every profile level  
     585            !-------------------------------------------------------------  
     586              
     587            DO ikn=1,inum_obs  
     588               iend=ista+ikn-1 
     589                   
     590               zweig(:,:,1) = 0._wp  
     591    
     592               ! This code forces the horizontal weights to be   
     593               ! zero IF the observation is below the bottom of the   
     594               ! corners of the interpolation nodes, Or if it is in   
     595               ! the mask. This is important for observations near   
     596               ! steep bathymetry  
     597               DO iin=1,2  
     598                  DO ijn=1,2  
     599      
     600                     depth_loop2: DO ik=kpk,2,-1  
     601                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     602                             
     603                           zweig(iin,ijn,1) = &   
     604                              & zweig2(iin,ijn,1) * &  
     605                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     606                              &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     607                             
     608                           EXIT depth_loop2  
     609 
     610                        ENDIF  
     611 
     612                     ENDDO depth_loop2  
     613      
     614                  ENDDO  
     615               ENDDO  
     616    
     617               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     618                  &              prodatqc%var(2)%vmod(iend:iend) )  
     619 
     620                  ! Set QC flag for any observations found below the bottom 
     621                  ! needed as the check here is more strict than that in obs_prep 
     622               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     623  
     624            ENDDO  
     625  
     626            DEALLOCATE(interp_corner,iv_indic)  
     627           
     628         ENDIF  
    447629 
    448630      ! Deallocate the data for interpolation 
    … …  
    459641         & zmask2, & 
    460642         & zint1,  & 
    461          & zint2   & 
     643         & zint2,  & 
     644         & zgdept, & 
     645         & zgdepw  & 
    462646         & ) 
    463647 

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

    • Property svn:keywords deleted
    r5785 r7773  
    2424   USE obs_inter_sup      ! Interpolation support 
    2525   USE obs_oper           ! Observation operators 
     26#if defined key_bdy 
     27   USE bdy_oce, ONLY : &        ! Boundary information 
     28      idx_bdy, nb_bdy 
     29#endif 
    2630   USE lib_mpp, ONLY : & 
    2731      & ctl_warn, ctl_stop 
    … …  
    4549CONTAINS 
    4650 
    47    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     51   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject ) 
    4852      !!---------------------------------------------------------------------- 
    4953      !!                    ***  ROUTINE obs_pre_sla  *** 
    … …  
    7276      !! * Arguments 
    7377      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data 
    74       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    75       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
     78      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening 
     79      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
     80      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
    7681      !! * Local declarations 
    7782      INTEGER :: iyea0        ! Initial date 
    … …  
    8792      INTEGER :: inlasobs     !  - close to land 
    8893      INTEGER :: igrdobs      !  - fail the grid search 
     94      INTEGER :: ibdysobs     !  - close to open boundary 
    8995                              ! Global counters for observations that 
    9096      INTEGER :: iotdobsmpp     !  - outside time domain 
    … …  
    9399      INTEGER :: inlasobsmpp    !  - close to land 
    94100      INTEGER :: igrdobsmpp     !  - fail the grid search 
     101      INTEGER :: ibdysobsmpp  !  - close to open boundary 
    95102      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    96103         & llvalid            ! SLA data selection 
    … …  
    118125      ilansobs = 0 
    119126      inlasobs = 0 
     127      ibdysobs = 0  
    120128 
    121129      ! ----------------------------------------------------------------------- 
    … …  
    151159         &                 tmask(:,:,1), surfdata%nqc,  & 
    152160         &                 iosdsobs,     ilansobs,     & 
    153          &                 inlasobs,     ld_nea        ) 
     161         &                 inlasobs,     ld_nea,       & 
     162         &                 ibdysobs,     ld_bound_reject        ) 
    154163 
    155164      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    156165      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    157166      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     167      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    158168 
    159169      ! ----------------------------------------------------------------------- 
    … …  
    201211               &            inlasobsmpp 
    202212         ENDIF 
     213         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     214            &            ibdysobsmpp   
    203215         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    204216            &            surfdataqc%nsurfmpp 
    … …  
    236248      &                     kpi, kpj, kpk, & 
    237249      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    238       &                     ld_nea, kdailyavtypes ) 
     250      &                     ld_nea, ld_bound_reject, kdailyavtypes ) 
    239251 
    240252!!---------------------------------------------------------------------- 
    … …  
    265277      LOGICAL, INTENT(IN) :: ld_var2 
    266278      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     279      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    267280      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    268281      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    … …  
    292305      INTEGER :: inlav1obs    !  - close to land (variable 1) 
    293306      INTEGER :: inlav2obs    !  - close to land (variable 2) 
     307      INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
     308      INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    294309      INTEGER :: igrdobs      !  - fail the grid search 
    295310      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
    … …  
    303318      INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    304319      INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
     320      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
     321      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    305322      INTEGER :: igrdobsmpp   !  - fail the grid search 
    306323      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
    … …  
    328345      ! Diagnotics counters for various failures. 
    329346 
    330       iotdobs  = 0 
    331       igrdobs  = 0 
     347      iotdobs   = 0 
     348      igrdobs   = 0 
    332349      iosdv1obs = 0 
    333350      iosdv2obs = 0 
    … …  
    336353      inlav1obs = 0 
    337354      inlav2obs = 0 
    338       iuvchku  = 0 
    339       iuvchkv = 0 
     355      ibdyv1obs = 0 
     356      ibdyv2obs = 0 
     357      iuvchku   = 0 
     358      iuvchkv   = 0 
    340359 
    341360      ! ----------------------------------------------------------------------- 
    … …  
    395414         &                 gdept_1d,              zmask1,               & 
    396415         &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    397          &                 iosdv1obs,              ilanv1obs,           & 
    398          &                 inlav1obs,              ld_nea                ) 
     416         &                 iosdv1obs,             ilanv1obs,            & 
     417         &                 inlav1obs,             ld_nea,               & 
     418         &                 ibdyv1obs,             ld_bound_reject       ) 
    399419 
    400420      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    401421      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    402422      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
     423      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    403424 
    404425      ! Variable 2 
    … …  
    414435         &                 gdept_1d,              zmask2,               & 
    415436         &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    416          &                 iosdv2obs,              ilanv2obs,           & 
    417          &                 inlav2obs,              ld_nea                ) 
     437         &                 iosdv2obs,             ilanv2obs,            & 
     438         &                 inlav2obs,             ld_nea,               & 
     439         &                 ibdyv2obs,             ld_bound_reject       ) 
    418440 
    419441      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    420442      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    421443      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     444      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
    422445 
    423446      ! ----------------------------------------------------------------------- 
    … …  
    489512               &            iuvchku 
    490513         ENDIF 
     514         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
     515               &            ibdyv1obsmpp 
    491516         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    492517            &            prodatqc%nvprotmpp(1) 
    … …  
    506531               &            iuvchkv 
    507532         ENDIF 
     533         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     534               &            ibdyv2obsmpp 
    508535         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    509536            &            prodatqc%nvprotmpp(2) 
    … …  
    875902      &                       plam,   pphi,    pmask,            & 
    876903      &                       kobsqc, kosdobs, klanobs,          & 
    877       &                       knlaobs,ld_nea                     ) 
     904      &                       knlaobs,ld_nea,                    & 
     905      &                       kbdyobs,ld_bound_reject            ) 
    878906      !!---------------------------------------------------------------------- 
    879907      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
    … …  
    908936      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    909937         & kobsqc             ! Observation quality control 
    910       INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain 
    911       INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell 
    912       INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land 
    913       LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land 
     938      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     939      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     940      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     941      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     942      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     943      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
    914944      !! * Local declarations 
    915945      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    916946         & zgmsk              ! Grid mask 
     947#if defined key_bdy  
     948      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     949         & zbmsk              ! Boundary mask 
     950      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     951#endif  
    917952      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    918953         & zglam, &           ! Model longitude at grid points 
    … …  
    956991 
    957992      END DO 
     993 
     994#if defined key_bdy              
     995      ! Create a mask grid points in boundary rim 
     996      IF (ld_bound_reject) THEN 
     997         zbdymask(:,:) = 1.0_wp 
     998         DO ji = 1, nb_bdy 
     999            DO jj = 1, idx_bdy(ji)%nblen(1) 
     1000               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1001            ENDDO 
     1002         ENDDO 
     1003  
     1004         CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, zbdymask, zbmsk )        
     1005      ENDIF 
     1006#endif        
    9581007       
    9591008      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
    … …  
    10001049            END DO 
    10011050         END DO 
    1002    
    1003          ! For observations on the grid reject them if their are at 
    1004          ! a masked point 
    1005           
     1051  
    10061052         IF (lgridobs) THEN 
    10071053            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    … …  
    10111057            ENDIF 
    10121058         ENDIF 
    1013                        
     1059 
     1060  
    10141061         ! Flag if the observation falls is close to land 
    10151062         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1016             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10171063            knlaobs = knlaobs + 1 
    1018             CYCLE 
    1019          ENDIF 
     1064            IF (ld_nea) THEN 
     1065               kobsqc(jobs) = kobsqc(jobs) + 14 
     1066               CYCLE 
     1067            ENDIF 
     1068         ENDIF 
     1069 
     1070#if defined key_bdy 
     1071         ! Flag if the observation falls close to the boundary rim 
     1072         IF (ld_bound_reject) THEN 
     1073            IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1074               kobsqc(jobs) = kobsqc(jobs) + 15 
     1075               kbdyobs = kbdyobs + 1 
     1076               CYCLE 
     1077            ENDIF 
     1078            ! for observations on the grid... 
     1079            IF (lgridobs) THEN 
     1080               IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1081                  kobsqc(jobs) = kobsqc(jobs) + 15 
     1082                  kbdyobs = kbdyobs + 1 
     1083                  CYCLE 
     1084               ENDIF 
     1085            ENDIF 
     1086         ENDIF 
     1087#endif  
    10201088             
    10211089      END DO 
    … …  
    10291097      &                       plam,    pphi,    pdep,    pmask, & 
    10301098      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1031       &                       klanobs, knlaobs, ld_nea          ) 
     1099      &                       klanobs, knlaobs, ld_nea,         & 
     1100      &                       kbdyobs, ld_bound_reject          ) 
    10321101      !!---------------------------------------------------------------------- 
    10331102      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
    … …  
    10521121      !! * Modules used 
    10531122      USE dom_oce, ONLY : &       ! Geographical information 
    1054          & gdepw_1d                         
     1123         & gdepw_1d,      & 
     1124         & gdepw_0,       &                        
     1125#if defined key_vvl 
     1126         & gdepw_n,       &  
     1127         & gdept_n,       & 
     1128#endif 
     1129         & ln_zco,        & 
     1130         & ln_zps,        & 
     1131         & lk_vvl                         
    10551132 
    10561133      !! * Arguments 
    … …  
    10861163      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10871164      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1165      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10881166      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1167      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
    10891168      !! * Local declarations 
    10901169      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10911170         & zgmsk              ! Grid mask 
     1171#if defined key_bdy  
     1172      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1173         & zbmsk              ! Boundary mask 
     1174      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     1175#endif  
     1176      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
     1177         & zgdepw          
    10921178      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
    10931179         & zglam, &           ! Model longitude at grid points 
    … …  
    10971183         & igrdj 
    10981184      LOGICAL :: lgridobs           ! Is observation on a model grid point. 
     1185      LOGICAL :: ll_next_to_land    ! Is a profile next to land  
    10991186      INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    11001187      INTEGER :: jobs, jobsp, jk, ji, jj 
    … …  
    11311218          
    11321219      END DO 
     1220 
     1221#if defined key_bdy  
     1222      ! Create a mask grid points in boundary rim 
     1223      IF (ld_bound_reject) THEN            
     1224         zbdymask(:,:) = 1.0_wp 
     1225         DO ji = 1, nb_bdy 
     1226            DO jj = 1, idx_bdy(ji)%nblen(1) 
     1227               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1228            ENDDO 
     1229         ENDDO 
     1230      ENDIF 
     1231  
     1232      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, zbdymask, zbmsk ) 
     1233#endif  
    11331234       
    11341235      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    … …  
    11591260         END DO 
    11601261 
     1262         ! Check if next to land  
     1263         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN  
     1264            ll_next_to_land=.TRUE.  
     1265         ELSE  
     1266            ll_next_to_land=.FALSE.  
     1267         ENDIF  
     1268          
    11611269         ! Reject observations 
    11621270 
    … …  
    11751283            ENDIF 
    11761284 
    1177             ! Flag if the observation falls with a model land cell 
    1178             IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1179                &  == 0.0_wp ) THEN 
    1180                kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1181                klanobs = klanobs + 1 
    1182                CYCLE 
     1285            ! To check if an observations falls within land there are two cases:  
     1286            ! 1: z-coordibnates, where the check uses the mask  
     1287            ! 2: terrain following (eg s-coordinates),   
     1288            !    where we use the depth of the bottom cell to mask observations  
     1289              
     1290            IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1)  
     1291                 
     1292               ! Flag if the observation falls with a model land cell  
     1293               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &  
     1294                  &  == 0.0_wp ) THEN  
     1295                  kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1296                  klanobs = klanobs + 1  
     1297                  CYCLE  
     1298               ENDIF  
     1299              
     1300               ! Flag if the observation is close to land  
     1301               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &  
     1302                  &  0.0_wp) THEN  
     1303                  knlaobs = knlaobs + 1  
     1304                  IF (ld_nea) THEN    
     1305                     kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1306                  ENDIF   
     1307               ENDIF  
     1308              
     1309            ELSE ! Case 2  
     1310               ! Flag if the observation is deeper than the bathymetry  
     1311               ! Or if it is within the mask  
     1312               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
     1313                  &     .OR. &  
     1314                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
     1315                  &  == 0.0_wp) ) THEN 
     1316                  kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1317                  klanobs = klanobs + 1  
     1318                  CYCLE  
     1319               ENDIF  
     1320                 
     1321               ! Flag if the observation is close to land  
     1322               IF ( ll_next_to_land ) THEN  
     1323                  knlaobs = knlaobs + 1  
     1324                  IF (ld_nea) THEN    
     1325                     kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1326                  ENDIF   
     1327               ENDIF  
     1328              
    11831329            ENDIF 
    11841330 
    … …  
    11941340            ENDIF 
    11951341             
    1196             ! Flag if the observation falls is close to land 
    1197             IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
    1198                &  0.0_wp) THEN 
    1199                IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14 
    1200                knlaobs = knlaobs + 1 
    1201             ENDIF 
    1202  
    12031342            ! Set observation depth equal to that of the first model depth 
    12041343            IF ( pobsdep(jobsp) <= pdep(1) ) THEN 
    12051344               pobsdep(jobsp) = pdep(1)   
    12061345            ENDIF 
     1346             
     1347#if defined key_bdy 
     1348            ! Flag if the observation falls close to the boundary rim 
     1349            IF (ld_bound_reject) THEN 
     1350               IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1351                  kobsqc(jobsp) = kobsqc(jobsp) + 15 
     1352                  kbdyobs = kbdyobs + 1 
     1353                  CYCLE 
     1354               ENDIF 
     1355               ! for observations on the grid... 
     1356               IF (lgridobs) THEN 
     1357                  IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1358                     kobsqc(jobsp) = kobsqc(jobsp) + 15 
     1359                     kbdyobs = kbdyobs + 1 
     1360                     CYCLE 
     1361                  ENDIF 
     1362               ENDIF 
     1363            ENDIF 
     1364#endif  
    12071365             
    12081366         END DO 

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted
    r5704 r7773  
    4444   !!---------------------------------------------------------------------- 
    4545   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    46    !! $Id$ 
     46   !! $Id: obs_read_altbias.F90 5704 2015-08-21 13:00:38Z mattmartin $ 
    4747   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4848   !!---------------------------------------------------------------------- 
    … …  
    128128         ! Get the Alt bias data 
    129129          
    130          CALL iom_get( numaltbias, jpdom_data, 'altbias', z_altbias(:,:), 1 ) 
     130         CALL iom_get( numaltbias, jpdom_autoglo, 'altbias', z_altbias(:,:), 1 ) 
    131131          
    132132         ! Close the file 

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    r7740 r7773  
    11MODULE obs_sstbias 
    22   !!====================================================================== 
    3    !!                       ***  MODULE obs_readsstbias  *** 
    4    !! Observation diagnostics: Read the bias for SLA data 
     3   !!                       ***  MODULE obs_sstbias  *** 
     4   !! Observation diagnostics: Read the bias for SST data 
    55   !!====================================================================== 
    66   !!---------------------------------------------------------------------- 
    7    !!   obs_rea_sstbias : Driver for reading altimeter bias 
     7   !!   obs_app_sstbias : Driver for reading and applying the SST bias 
    88   !!---------------------------------------------------------------------- 
    99   !! * Modules used    
    … …  
    2222   USE dom_oce, ONLY : &        ! Domain variables 
    2323      & tmask, & 
    24       & tmask_i, & 
    25       & e1t,   & 
    26       & e2t,   & 
    2724      & gphit, & 
    2825      & glamt 
    29    USE oce, ONLY : &           ! Model variables 
    30       & sshn 
    3126   USE obs_inter_h2d 
    3227   USE obs_utils               ! Various observation tools 
    … …  
    3732   PUBLIC obs_app_sstbias     ! Read the altimeter bias 
    3833CONTAINS 
    39    SUBROUTINE obs_app_sstbias( ksstno, sstdata, k2dint, knumtypes, & 
     34   SUBROUTINE obs_app_sstbias( sstdata, k2dint, knumtypes, & 
    4035                               cl_bias_files ) 
    4136      !!--------------------------------------------------------------------- 
    4237      !! 
    43       !!                   *** ROUTINE obs_rea_sstbias *** 
     38      !!                   *** ROUTINE obs_app_sstbias *** 
    4439      !! 
    4540      !! ** Purpose : Read SST bias data from files and apply correction to 
    … …  
    5954      USE iom 
    6055      USE netcdf 
     56 
    6157      !! * Arguments 
    62       INTEGER, INTENT(IN) :: ksstno      ! Number of SST obs sets 
    63       TYPE(obs_surf), DIMENSION(ksstno), INTENT(INOUT) :: & 
    64          & sstdata       ! SST data 
     58      TYPE(obs_surf), INTENT(INOUT) :: & 
     59         & sstdata            ! SST data 
    6560      INTEGER, INTENT(IN) :: k2dint 
    66       INTEGER, INTENT(IN) :: knumtypes !number of bias types to read in 
     61      INTEGER, INTENT(IN) :: & 
     62         & knumtypes          ! Number of bias types to read in 
    6763      CHARACTER(LEN=128), DIMENSION(knumtypes), INTENT(IN) :: & 
    68                           cl_bias_files !List of files to read 
     64         & cl_bias_files      ! List of files to read 
     65 
    6966      !! * Local declarations 
    7067      INTEGER :: jslano       ! Data set loop variable 
    … …  
    8077      INTEGER :: i_var_id 
    8178      INTEGER, DIMENSION(knumtypes) :: & 
    82          & ibiastypes             ! Array of the bias types in each file 
     79         & ibiastypes         ! Array of the bias types in each file 
    8380      REAL(wp), DIMENSION(jpi,jpj,knumtypes) :: &  
    84          & z_sstbias              ! Array to store the SST bias values 
     81         & z_sstbias          ! Array to store the SST bias values 
    8582      REAL(wp), DIMENSION(jpi,jpj) :: &  
    86          & z_sstbias_2d           ! Array to store the SST bias values    
     83         & z_sstbias_2d       ! Array to store the SST bias values    
    8784      REAL(wp), DIMENSION(1) :: & 
    8885         & zext, & 
    … …  
    114111      INTEGER :: iret  
    115112      INTEGER :: inumtype 
    116       IF(lwp)WRITE(numout,*)  
    117       IF(lwp)WRITE(numout,*) 'obs_rea_sstbias : ' 
    118       IF(lwp)WRITE(numout,*) '----------------- ' 
    119       IF(lwp)WRITE(numout,*) 'Read SST bias ' 
    120       ! Open and read the files 
    121       z_sstbias(:,:,:)=0.0_wp 
     113 
     114      IF ( lwp ) THEN 
     115         WRITE(numout,*)  
     116         WRITE(numout,*) 'obs_app_sstbias : ' 
     117         WRITE(numout,*) '----------------- ' 
     118         WRITE(numout,*) 'Read SST bias ' 
     119      ENDIF 
     120 
     121      ! Open and read the SST bias files for each bias type 
     122      z_sstbias(:,:,:) = 0.0_wp 
     123 
    122124      DO jtype = 1, knumtypes 
    123125      
    124126         numsstbias=0 
    125          IF(lwp)WRITE(numout,*) 'Opening ',cl_bias_files(jtype) 
    126          CALL iom_open( cl_bias_files(jtype), numsstbias, ldstop=.FALSE. )        
     127 
     128         IF ( lwp ) WRITE(numout,*) 'Opening ',cl_bias_files(jtype) 
     129         CALL iom_open( cl_bias_files(jtype), numsstbias, ldstop=.FALSE. ) 
     130 
    127131         IF (numsstbias .GT. 0) THEN 
    128132      
    … …  
    137141            iret=NF90_CLOSE(incfile)        
    138142            
    139             IF ( iret /= 0  ) CALL ctl_stop( & 
    140                'obs_rea_sstbias : Cannot read bias type from file '// & 
    141                cl_bias_files(jtype) ) 
     143            IF ( iret /= 0  ) THEN 
     144               CALL ctl_stop( 'obs_app_sstbias : Cannot read bias type from file '// & 
     145                  &           TRIM( cl_bias_files(jtype) ) ) 
     146            ENDIF 
     147 
    142148            ! Get the SST bias data 
    143149            CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 ) 
    144150            z_sstbias(:,:,jtype) = z_sstbias_2d(:,:)        
    145151            ! Close the file 
    146             CALL iom_close(numsstbias)        
     152            CALL iom_close(numsstbias) 
     153      
    147154         ELSE 
    148155            CALL ctl_stop('obs_read_sstbias: File '// &  
    149                            TRIM( cl_bias_files(jtype) )//' Not found') 
     156               &          TRIM( cl_bias_files(jtype) )//' Not found') 
    150157         ENDIF 
     158 
    151159      END DO 
    152160            
    153       ! Interpolate the bias already on the model grid at the observation point 
    154       DO jslano = 1, ksstno 
     161      ! Interpolate the bias from the model grid to the observation points 
     162      ALLOCATE( & 
     163         & igrdi(2,2,sstdata%nsurf), & 
     164         & igrdj(2,2,sstdata%nsurf), & 
     165         & zglam(2,2,sstdata%nsurf), & 
     166         & zgphi(2,2,sstdata%nsurf), & 
     167         & zmask(2,2,sstdata%nsurf)  ) 
     168        
     169      DO jobs = 1, sstdata%nsurf  
     170         igrdi(1,1,jobs) = sstdata%mi(jobs)-1 
     171         igrdj(1,1,jobs) = sstdata%mj(jobs)-1 
     172         igrdi(1,2,jobs) = sstdata%mi(jobs)-1 
     173         igrdj(1,2,jobs) = sstdata%mj(jobs) 
     174         igrdi(2,1,jobs) = sstdata%mi(jobs) 
     175         igrdj(2,1,jobs) = sstdata%mj(jobs)-1 
     176         igrdi(2,2,jobs) = sstdata%mi(jobs) 
     177         igrdj(2,2,jobs) = sstdata%mj(jobs) 
     178      END DO 
     179 
     180      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, & 
     181         &                  igrdi, igrdj, glamt, zglam ) 
     182      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, & 
     183         &                  igrdi, igrdj, gphit, zgphi ) 
     184      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, & 
     185         &                  igrdi, igrdj, tmask(:,:,1), zmask ) 
     186 
     187      DO jtype = 1, knumtypes 
     188          
     189         !Find the number observations of type 
     190         !and alllocate tempory arrays 
     191         inumtype = COUNT( sstdata%ntyp(:) == ibiastypes(jtype) ) 
     192 
    155193         ALLOCATE( & 
    156             & igrdi(2,2,sstdata(jslano)%nsurf), & 
    157             & igrdj(2,2,sstdata(jslano)%nsurf), & 
    158             & zglam(2,2,sstdata(jslano)%nsurf), & 
    159             & zgphi(2,2,sstdata(jslano)%nsurf), & 
    160             & zmask(2,2,sstdata(jslano)%nsurf)  ) 
    161         
    162          DO jobs = 1, sstdata(jslano)%nsurf  
    163             igrdi(1,1,jobs) = sstdata(jslano)%mi(jobs)-1 
    164             igrdj(1,1,jobs) = sstdata(jslano)%mj(jobs)-1 
    165             igrdi(1,2,jobs) = sstdata(jslano)%mi(jobs)-1 
    166             igrdj(1,2,jobs) = sstdata(jslano)%mj(jobs) 
    167             igrdi(2,1,jobs) = sstdata(jslano)%mi(jobs) 
    168             igrdj(2,1,jobs) = sstdata(jslano)%mj(jobs)-1 
    169             igrdi(2,2,jobs) = sstdata(jslano)%mi(jobs) 
    170             igrdj(2,2,jobs) = sstdata(jslano)%mj(jobs) 
    171          END DO 
    172          CALL obs_int_comm_2d( 2, 2, sstdata(jslano)%nsurf, & 
    173             &                  igrdi, igrdj, glamt, zglam ) 
    174          CALL obs_int_comm_2d( 2, 2, sstdata(jslano)%nsurf, & 
    175             &                  igrdi, igrdj, gphit, zgphi ) 
    176          CALL obs_int_comm_2d( 2, 2, sstdata(jslano)%nsurf, & 
    177             &                  igrdi, igrdj, tmask(:,:,1), zmask ) 
    178          DO jtype = 1, knumtypes 
    179           
    180             !Find the number observations of type 
    181             !and alllocate tempory arrays 
    182             inumtype = COUNT( sstdata(jslano)%ntyp(:) == ibiastypes(jtype) ) 
    183             ALLOCATE( & 
    184194            & igrdi_tmp(2,2,inumtype), & 
    185195            & igrdj_tmp(2,2,inumtype), & 
    … …  
    188198            & zmask_tmp(2,2,inumtype), & 
    189199            & zbias( 2,2,inumtype ) ) 
    190             jt=1 
    191             DO jobs = 1, sstdata(jslano)%nsurf  
    192                IF ( sstdata(jslano)%ntyp(jobs) == ibiastypes(jtype) ) THEN 
    193                   igrdi_tmp(:,:,jt) = igrdi(:,:,jobs)  
    194                   igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 
    195                   zglam_tmp(:,:,jt) = zglam(:,:,jobs) 
    196                   zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    197                   zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    198                   zmask_tmp(:,:,jt) = zmask(:,:,jobs) 
    199                   jt = jt +1 
    200                ENDIF 
    201             END DO 
     200 
     201         jt=1 
     202         DO jobs = 1, sstdata%nsurf  
     203            IF ( sstdata%ntyp(jobs) == ibiastypes(jtype) ) THEN 
     204 
     205               igrdi_tmp(:,:,jt) = igrdi(:,:,jobs)  
     206               igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 
     207               zglam_tmp(:,:,jt) = zglam(:,:,jobs) 
     208               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
     209               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
     210               zmask_tmp(:,:,jt) = zmask(:,:,jobs) 
     211 
     212               jt = jt +1 
     213 
     214            ENDIF 
     215         END DO 
    202216                          
    203             CALL obs_int_comm_2d( 2, 2, inumtype, & 
    204                   &           igrdi_tmp(:,:,:), igrdj_tmp(:,:,:), & 
    205                   &           z_sstbias(:,:,jtype), zbias(:,:,:) ) 
    206             jt=1 
    207             DO jobs = 1, sstdata(jslano)%nsurf 
    208                IF ( sstdata(jslano)%ntyp(jobs) == ibiastypes(jtype) ) THEN 
    209                   zlam = sstdata(jslano)%rlam(jobs) 
    210                   zphi = sstdata(jslano)%rphi(jobs) 
    211                   iico = sstdata(jslano)%mi(jobs) 
    212                   ijco = sstdata(jslano)%mj(jobs)          
    213                   CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    214                      &                   zglam_tmp(:,:,jt), & 
    215                      &                   zgphi_tmp(:,:,jt), & 
    216                      &                   zmask_tmp(:,:,jt), zweig, zobsmask ) 
    217                   CALL obs_int_h2d( 1, 1,      & 
    218                      &              zweig, zbias(:,:,jt),  zext ) 
    219                   ! adjust sst with bias field 
    220                   sstdata(jslano)%robs(jobs,1) = & 
    221                      sstdata(jslano)%robs(jobs,1) - zext(1) 
    222                   jt=jt+1 
    223                ENDIF 
    224             END DO  
     217         CALL obs_int_comm_2d( 2, 2, inumtype, & 
     218               &           igrdi_tmp(:,:,:), igrdj_tmp(:,:,:), & 
     219               &           z_sstbias(:,:,jtype), zbias(:,:,:) ) 
     220 
     221         jt=1 
     222         DO jobs = 1, sstdata%nsurf 
     223            IF ( sstdata%ntyp(jobs) == ibiastypes(jtype) ) THEN 
     224 
     225               zlam = sstdata%rlam(jobs) 
     226               zphi = sstdata%rphi(jobs) 
     227               iico = sstdata%mi(jobs) 
     228               ijco = sstdata%mj(jobs)          
     229 
     230               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     231                  &                   zglam_tmp(:,:,jt), & 
     232                  &                   zgphi_tmp(:,:,jt), & 
     233                  &                   zmask_tmp(:,:,jt), zweig, zobsmask ) 
     234 
     235               CALL obs_int_h2d( 1, 1, zweig, zbias(:,:,jt),  zext ) 
     236 
     237               ! adjust sst with bias field 
     238               sstdata%robs(jobs,1) = & 
     239                  &    sstdata%robs(jobs,1) - zext(1) 
     240 
     241               jt=jt+1 
     242 
     243            ENDIF 
     244         END DO  
    225245                
    226             !Deallocate arrays 
    227             DEALLOCATE( & 
    228             & igrdi_tmp, & 
    229             & igrdj_tmp, & 
    230             & zglam_tmp, & 
    231             & zgphi_tmp, & 
    232             & zmask_tmp, & 
    233             & zbias )            
    234          END DO 
     246         !Deallocate arrays 
    235247         DEALLOCATE( & 
    236             & igrdi, & 
    237             & igrdj, & 
    238             & zglam, & 
    239             & zgphi, & 
    240             & zmask ) 
    241       END DO 
     248         & igrdi_tmp, & 
     249         & igrdj_tmp, & 
     250         & zglam_tmp, & 
     251         & zgphi_tmp, & 
     252         & zmask_tmp, & 
     253         & zbias )       
     254      
     255      END DO !jtype 
     256 
     257      DEALLOCATE( & 
     258         & igrdi, & 
     259         & igrdj, & 
     260         & zglam, & 
     261         & zgphi, & 
     262         & zmask ) 
     263 
    242264      IF(lwp) THEN 
    243265         WRITE(numout,*) " " 
    244266         WRITE(numout,*) "SST bias correction applied successfully" 
    245267         WRITE(numout,*) "Obs types: ",ibiastypes(:), & 
    246                               " Have all been bias corrected\n" 
     268                              " have all been bias corrected\n" 
    247269      ENDIF 
    248270   END SUBROUTINE obs_app_sstbias 

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

    • Property svn:keywords deleted
    r5682 r7773  
    5050      INTEGER :: npj 
    5151      INTEGER :: nsurfup    !: Observation counter used in obs_oper 
     52      INTEGER :: nrec       !: Number of surface observation records in window 
    5253 
    5354      ! Arrays with size equal to the number of surface observations 
    … …  
    5657         & mi,   &        !: i-th grid coord. for interpolating to surface observation 
    5758         & mj,   &        !: j-th grid coord. for interpolating to surface observation 
     59         & mt,   &        !: time record number for gridded data 
    5860         & nsidx,&        !: Surface observation number 
    5961         & nsfil,&        !: Surface observation number in file 
    … …  
    9395         & nsstpmpp       !: Global number of surface observations per time step 
    9496 
     97      ! Arrays with size equal to the number of observation records in the window 
     98      INTEGER, POINTER, DIMENSION(:) :: & 
     99         & mrecstp   ! Time step of the records 
     100 
    95101      ! Arrays used to store source indices when  
    96102      ! compressing obs_surf derived types 
    … …  
    101107         & nsind          !: Source indices of surface data in compressed data 
    102108 
     109      ! Is this a gridded product? 
     110      
     111      LOGICAL :: lgrid 
     112 
    103113   END TYPE obs_surf 
    104114 
    105115   !!---------------------------------------------------------------------- 
    106116   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    107    !! $Id$ 
     117   !! $Id: obs_surf_def.F90 5682 2015-08-12 15:46:45Z mattmartin $ 
    108118   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    109119   !!---------------------------------------------------------------------- 
    … …  
    160170         & surf%mi(ksurf),      & 
    161171         & surf%mj(ksurf),      & 
     172         & surf%mt(ksurf),      & 
    162173         & surf%nsidx(ksurf),   & 
    163174         & surf%nsfil(ksurf),   & 
    … …  
    176187         & ) 
    177188 
     189      surf%mt(:) = -1 
     190 
    178191 
    179192      ! Allocate arrays of number of surface data size * number of variables 
    … …  
    190203         & ) 
    191204 
     205      surf%rext(:,:) = 0.0_wp  
     206 
    192207      ! Allocate arrays of number of time step size 
    193208 
    … …  
    217232 
    218233      surf%nsurfup     = 0 
     234       
     235      ! Not gridded by default 
     236           
     237      surf%lgrid       = .FALSE. 
    219238               
    220239   END SUBROUTINE obs_surf_alloc 
    … …  
    242261         & surf%mi,      & 
    243262         & surf%mj,      & 
     263         & surf%mt,      & 
    244264         & surf%nsidx,   & 
    245265         & surf%nsfil,   & 
    … …  
    370390            newsurf%mi(insurf)    = surf%mi(ji) 
    371391            newsurf%mj(insurf)    = surf%mj(ji) 
     392            newsurf%mt(insurf)    = surf%mt(ji) 
    372393            newsurf%nsidx(insurf) = surf%nsidx(ji) 
    373394            newsurf%nsfil(insurf) = surf%nsfil(ji) 
    … …  
    414435      newsurf%nstp     = surf%nstp 
    415436      newsurf%cvars(:) = surf%cvars(:) 
     437       
     438      ! Set gridded stuff 
     439       
     440      newsurf%mt(insurf)    = surf%mt(ji) 
    416441  
    417442      ! Deallocate temporary data 
    … …  
    454479         oldsurf%mi(jj)    = surf%mi(ji) 
    455480         oldsurf%mj(jj)    = surf%mj(ji) 
     481         oldsurf%mt(jj)    = surf%mt(ji) 
    456482         oldsurf%nsidx(jj) = surf%nsidx(ji) 
    457483         oldsurf%nsfil(jj) = surf%nsfil(ji) 

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted

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

    • Property svn:keywords deleted
    r5704 r7773  
    88   !!   obs_wri_prof   : Write profile observations in feedback format 
    99   !!   obs_wri_surf   : Write surface observations in feedback format 
    10    !!   obs_wri_stats : Print basic statistics on the data being written out 
     10   !!   obs_wri_stats  : Print basic statistics on the data being written out 
    1111   !!---------------------------------------------------------------------- 
    1212 
    … …  
    4848   !!---------------------------------------------------------------------- 
    4949   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    50    !! $Id$ 
     50   !! $Id: obs_write.F90 5704 2015-08-21 13:00:38Z mattmartin $ 
    5151   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5252   !!---------------------------------------------------------------------- 
    … …  
    411411         fbdata%caddlong(1,1) = 'Model interpolated ICE' 
    412412         fbdata%caddunit(1,1) = 'Fraction' 
     413         fbdata%cgrid(1)      = 'T' 
     414         DO ja = 1, iadd 
     415            fbdata%caddname(1+ja) = padd%cdname(ja) 
     416            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     417            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     418         END DO 
     419 
     420      CASE('SSS') 
     421 
     422         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     423            &                 1 + iadd, iext, .TRUE. ) 
     424 
     425         clfiletype = 'sssfb' 
     426         fbdata%cname(1)      = surfdata%cvars(1) 
     427         fbdata%coblong(1)    = 'Sea surface salinity' 
     428         fbdata%cobunit(1)    = 'psu' 
     429         DO je = 1, iext 
     430            fbdata%cextname(je) = pext%cdname(je) 
     431            fbdata%cextlong(je) = pext%cdlong(je,1) 
     432            fbdata%cextunit(je) = pext%cdunit(je,1) 
     433         END DO 
     434         fbdata%caddlong(1,1) = 'Model interpolated SSS' 
     435         fbdata%caddunit(1,1) = 'psu' 
     436         fbdata%cgrid(1)      = 'T' 
     437         DO ja = 1, iadd 
     438            fbdata%caddname(1+ja) = padd%cdname(ja) 
     439            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     440            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     441         END DO 
     442 
     443      CASE('LOGCHL') 
     444 
     445         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     446            &                 1 + iadd, iext, .TRUE. ) 
     447 
     448         clfiletype = 'logchlfb' 
     449         fbdata%cname(1)      = surfdata%cvars(1) 
     450         fbdata%coblong(1)    = 'logchl concentration' 
     451         fbdata%cobunit(1)    = 'mg/m3' 
     452         DO je = 1, iext 
     453            fbdata%cextname(je) = pext%cdname(je) 
     454            fbdata%cextlong(je) = pext%cdlong(je,1) 
     455            fbdata%cextunit(je) = pext%cdunit(je,1) 
     456         END DO 
     457         fbdata%caddlong(1,1) = 'Model interpolated LOGCHL' 
     458         fbdata%caddunit(1,1) = 'mg/m3' 
     459         fbdata%cgrid(1)      = 'T' 
     460         DO ja = 1, iadd 
     461            fbdata%caddname(1+ja) = padd%cdname(ja) 
     462            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     463            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     464         END DO 
     465 
     466      CASE('SPM') 
     467 
     468         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     469            &                 1 + iadd, iext, .TRUE. ) 
     470 
     471         clfiletype = 'spmfb' 
     472         fbdata%cname(1)      = surfdata%cvars(1) 
     473         fbdata%coblong(1)    = 'spm' 
     474         fbdata%cobunit(1)    = 'g/m3' 
     475         DO je = 1, iext 
     476            fbdata%cextname(je) = pext%cdname(je) 
     477            fbdata%cextlong(je) = pext%cdlong(je,1) 
     478            fbdata%cextunit(je) = pext%cdunit(je,1) 
     479         END DO 
     480         fbdata%caddlong(1,1) = 'Model interpolated spm' 
     481         fbdata%caddunit(1,1) = 'g/m3' 
     482         fbdata%cgrid(1)      = 'T' 
     483         DO ja = 1, iadd 
     484            fbdata%caddname(1+ja) = padd%cdname(ja) 
     485            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     486            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     487         END DO 
     488 
     489      CASE('FCO2') 
     490 
     491         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     492            &                 1 + iadd, iext, .TRUE. ) 
     493 
     494         clfiletype = 'fco2fb' 
     495         fbdata%cname(1)      = surfdata%cvars(1) 
     496         fbdata%coblong(1)    = 'fco2' 
     497         fbdata%cobunit(1)    = 'uatm' 
     498         DO je = 1, iext 
     499            fbdata%cextname(je) = pext%cdname(je) 
     500            fbdata%cextlong(je) = pext%cdlong(je,1) 
     501            fbdata%cextunit(je) = pext%cdunit(je,1) 
     502         END DO 
     503         fbdata%caddlong(1,1) = 'Model interpolated fco2' 
     504         fbdata%caddunit(1,1) = 'uatm' 
     505         fbdata%cgrid(1)      = 'T' 
     506         DO ja = 1, iadd 
     507            fbdata%caddname(1+ja) = padd%cdname(ja) 
     508            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     509            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     510         END DO 
     511 
     512      CASE('PCO2') 
     513 
     514         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     515            &                 1 + iadd, iext, .TRUE. ) 
     516 
     517         clfiletype = 'pco2fb' 
     518         fbdata%cname(1)      = surfdata%cvars(1) 
     519         fbdata%coblong(1)    = 'pco2' 
     520         fbdata%cobunit(1)    = 'uatm' 
     521         DO je = 1, iext 
     522            fbdata%cextname(je) = pext%cdname(je) 
     523            fbdata%cextlong(je) = pext%cdlong(je,1) 
     524            fbdata%cextunit(je) = pext%cdunit(je,1) 
     525         END DO 
     526         fbdata%caddlong(1,1) = 'Model interpolated pco2' 
     527         fbdata%caddunit(1,1) = 'uatm' 
    413528         fbdata%cgrid(1)      = 'T' 
    414529         DO ja = 1, iadd 

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90

    • Property svn:keywords deleted
    r2474 r7773  
    11   !!---------------------------------------------------------------------- 
    22   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3    !! $Id$ 
     3   !! $Id: obsinter_h2d.h90 2474 2010-12-16 15:32:33Z djlea $ 
    44   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    55   !!---------------------------------------------------------------------- 
    … …  
    12401240         & zdum,  & 
    12411241         & zaamax 
    1242         
     1242       
     1243      imax = -1  
    12431244      ! Main computation 
    12441245      pflt = 1.0_wp 

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

    • Property svn:keywords deleted

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/str_c_to_for.h90

    • Property svn:keywords deleted
Note: See TracChangeset for help on using the changeset viewer.