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 15764 for branches/UKMO/r6232_obs_oper_update – NEMO

Ignore:
Timestamp:
2022-03-21T13:12:29+01:00 (2 years ago)
Author:
jcastill
Message:

Update with the latest changes in branches/UKMO/dev_r5518_obs_oper_update@15400

Location:
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r11203 r15764  
    2121   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments 
    2222   !!   ssh_asm_inc    : Apply the SSH increment 
    23    !!   seaice_asm_inc : Apply the seaice increment 
     23   !!   seaice_asm_inc : Apply the sea ice concentration increment 
     24   !!   sit_asm_inc    : Apply the sea ice thickness increment 
    2425   !!---------------------------------------------------------------------- 
    2526   USE wrk_nemo         ! Memory Allocation 
    … …  
    4950   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments 
    5051   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
    51    PUBLIC   seaice_asm_inc !: Apply the seaice increment 
     52   PUBLIC   seaice_asm_inc !: Apply the seaice concentration increment 
     53   PUBLIC   sit_asm_inc    !: Apply the seaice thickness increment 
     54   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments 
    5255 
    5356#if defined key_asminc 
    … …  
    5760#endif 
    5861   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
     62   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields 
    5963   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6064   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
    … …  
    6266   LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments 
    6367   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    64    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment 
     68   LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE.   !: No sea ice concentration increment 
     69   LOGICAL, PUBLIC :: ln_sitinc = .FALSE.      !: No sea ice thickness increment 
     70   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments 
    6571   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    6672   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
    … …  
    8591   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    8692   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     93   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   sit_bkginc            ! Increment to the background sea ice thickness 
    8794 
    8895   !! * Substitutions 
    … …  
    127134      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    128135      !! 
    129       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     136      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg, ln_balwri,                & 
    130137         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    131138         &                 ln_asmdin, ln_asmiau,                           & 
    … …  
    138145      !----------------------------------------------------------------------- 
    139146      ln_seaiceinc = .FALSE. 
     147      ln_sitinc = .FALSE. 
    140148      ln_temnofreeze = .FALSE. 
    141149 
    … …  
    154162         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :' 
    155163         WRITE(numout,*) '~~~~~~~~~~~~' 
    156          WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
     164         WRITE(numout,*) '   Namelist nam_asminc : set assimilation increment parameters' 
    157165         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     166         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg 
     167         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri 
    158168         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    159169         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
    160170         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc 
     171         WRITE(numout,*) '      Logical switch for applying SIC increments               ln_seaiceinc = ', ln_seaiceinc 
     172         WRITE(numout,*) '      Logical switch for applying SIT increments               ln_sitinc = ', ln_sitinc 
    161173         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
    162          WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc 
    163174         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    164175         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
    … …  
    216227 
    217228      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    218            .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 
    219          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 
     229         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 
     230         &        ( ln_sitinc ).OR.( lk_bgcinc ) )) & 
     231         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
     232         &                ' ln_sitinc and ln_(bgc-variable)inc is set to .true.', & 
    220233         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    221234         &                ' Inconsistent options') 
    … …  
    226239 
    227240      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    228          &                     )  & 
    229          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 
     241         & .AND.( .NOT. ln_sitinc ).AND.( .NOT. lk_bgcinc ) )  & 
     242         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
     243         &                ' ln_sitinc and ln_(bgc-variable)inc are set to .false. :', & 
    230244         &                ' The assimilation increments are not applied') 
    231245 
    … …  
    325339      !-------------------------------------------------------------------- 
    326340 
    327       ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
    328       ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
    329       ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
    330       ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
    331       ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
    332       ALLOCATE( seaice_bkginc(jpi,jpj)) 
     341      IF ( ln_trainc ) THEN 
     342         ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
     343         ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
     344         t_bkginc(:,:,:) = 0.0 
     345         s_bkginc(:,:,:) = 0.0 
     346      ENDIF 
     347      IF ( ln_dyninc ) THEN  
     348         ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
     349         ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
     350         u_bkginc(:,:,:) = 0.0 
     351         v_bkginc(:,:,:) = 0.0 
     352      ENDIF 
     353      IF ( ln_sshinc ) THEN 
     354         ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
     355         ssh_bkginc(:,:) = 0.0 
     356      ENDIF 
     357      IF ( ln_seaiceinc ) THEN  
     358         ALLOCATE( seaice_bkginc(jpi,jpj)) 
     359         seaice_bkginc(:,:) = 0.0 
     360      ENDIF 
     361      IF ( ln_sitinc ) THEN  
     362         ALLOCATE( sit_bkginc(jpi,jpj)) 
     363         sit_bkginc(:,:) = 0.0 
     364      ENDIF 
    333365#if defined key_asminc 
    334366      ALLOCATE( ssh_iau(jpi,jpj)      ) 
    335 #endif 
    336       t_bkginc(:,:,:) = 0.0 
    337       s_bkginc(:,:,:) = 0.0 
    338       u_bkginc(:,:,:) = 0.0 
    339       v_bkginc(:,:,:) = 0.0 
    340       ssh_bkginc(:,:) = 0.0 
    341       seaice_bkginc(:,:) = 0.0 
    342 #if defined key_asminc 
    343367      ssh_iau(:,:)    = 0.0 
    344368#endif 
    345       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 
     369      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 
     370         &  .OR.( ln_sitinc ).OR.( lk_bgcinc ) ) THEN 
    346371 
    347372         !-------------------------------------------------------------------- 
    … …  
    406431            ! to allow for differences in masks 
    407432            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0 
     433         ENDIF  
     434 
     435         IF ( ln_sitinc ) THEN  
     436            CALL iom_get( inum, jpdom_autoglo, 'bckinsit', sit_bkginc, 1 )  
     437            ! Apply the masks  
     438            sit_bkginc(:,:) = sit_bkginc(:,:) * tmask(:,:,1)  
     439            ! Set missing increments to 0.0 rather than 1e+20  
     440            ! to allow for differences in masks  
     441            WHERE( ABS( sit_bkginc(:,:) ) > 1.0e+10 ) sit_bkginc(:,:) = 0.0  
    408442         ENDIF 
    409443 
    … …  
    766800      ! Perhaps the following call should be in step 
    767801      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment 
     802      IF   ( ln_sitinc  )      CALL sit_asm_inc ( kt )      ! apply sea ice thickness increment 
    768803      ! 
    769804   END SUBROUTINE tra_asm_inc 
    … …  
    923958   END SUBROUTINE ssh_asm_inc 
    924959 
     960   SUBROUTINE sit_asm_inc( kt, kindic ) 
     961      !!---------------------------------------------------------------------- 
     962      !!                    ***  ROUTINE sit_asm_inc  *** 
     963      !!           
     964      !! ** Purpose : Apply the sea ice thickness assimilation increment. 
     965      !! 
     966      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
     967      !! 
     968      !! ** Action  :  
     969      !! 
     970      !!---------------------------------------------------------------------- 
     971      IMPLICIT NONE 
     972      ! 
     973      INTEGER, INTENT(in)           ::   kt   ! Current time step 
     974      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
     975      ! 
     976      INTEGER  ::   it 
     977      REAL(wp) ::   zincwgt   ! IAU weight for current time step 
     978      !!---------------------------------------------------------------------- 
     979 
     980      IF ( ln_asmiau ) THEN 
     981 
     982         !-------------------------------------------------------------------- 
     983         ! Incremental Analysis Updating 
     984         !-------------------------------------------------------------------- 
     985 
     986         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     987 
     988            it = kt - nit000 + 1 
     989            zincwgt = wgtiau(it)      ! IAU weight for the current time step  
     990            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
     991            ! EF: Actually CICE is expecting a tendency so is divided by rdt below 
     992 
     993            IF(lwp) THEN 
     994               WRITE(numout,*)  
     995               WRITE(numout,*) 'sit_asm_inc : sea ice thick IAU at time step = ', & 
     996                  &  kt,' with IAU weight = ', wgtiau(it) 
     997               WRITE(numout,*) '~~~~~~~~~~~~' 
     998            ENDIF 
     999 
     1000#if defined key_cice && defined key_asminc 
     1001            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 
     1002            ndsit_da(:,:) = sit_bkginc(:,:) * zincwgt / rdt 
     1003#endif 
     1004 
     1005            IF ( kt == nitiaufin_r ) THEN 
     1006               DEALLOCATE( sit_bkginc ) 
     1007            ENDIF 
     1008 
     1009         ELSE 
     1010 
     1011#if defined key_cice && defined key_asminc 
     1012            ! Sea-ice thickness : CICE case. Zero ice increment tendency into CICE 
     1013            ndsit_da(:,:) = 0.0_wp 
     1014#endif 
     1015 
     1016         ENDIF 
     1017 
     1018      ELSEIF ( ln_asmdin ) THEN 
     1019 
     1020         !-------------------------------------------------------------------- 
     1021         ! Direct Initialization 
     1022         !-------------------------------------------------------------------- 
     1023 
     1024         IF ( kt == nitdin_r ) THEN 
     1025 
     1026            neuler = 0                    ! Force Euler forward step 
     1027 
     1028#if defined key_cice && defined key_asminc 
     1029            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 
     1030           ndsit_da(:,:) = sit_bkginc(:,:) / rdt 
     1031#endif 
     1032           IF ( .NOT. PRESENT(kindic) ) THEN 
     1033              DEALLOCATE( sit_bkginc ) 
     1034           END IF 
     1035 
     1036         ELSE 
     1037 
     1038#if defined key_cice && defined key_asminc 
     1039            ! Sea-ice thicnkness : CICE case. Zero ice thickness increment tendency into CICE  
     1040            ndsit_da(:,:) = 0.0_wp 
     1041#endif 
     1042          
     1043         ENDIF 
     1044 
     1045      ENDIF 
     1046 
     1047   END SUBROUTINE sit_asm_inc 
    9251048 
    9261049   SUBROUTINE seaice_asm_inc( kt, kindic ) 
    … …  
    9281051      !!                    ***  ROUTINE seaice_asm_inc  *** 
    9291052      !!           
    930       !! ** Purpose : Apply the sea ice assimilation increment. 
     1053      !! ** Purpose : Apply the sea ice concentration assimilation increment. 
    9311054      !! 
    9321055      !! ** Method  : Direct initialization or Incremental Analysis Updating. 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r11203 r15764  
    4545   USE agrif_opa_interp ! agrif 
    4646#endif 
    47 #if defined key_asminc    
    48    USE asminc          ! Assimilation increment 
    49 #endif 
    5047 
    5148   IMPLICIT NONE 
    … …  
    459456                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    460457      ENDIF 
    461 #if defined key_asminc 
    462       !                                         ! Include the IAU weighted SSH increment 
    463       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    464          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    465       ENDIF 
    466 #endif 
    467       !                                   !* Fill boundary data arrays for AGRIF 
    468       !                                   ! ------------------------------------ 
     458      !                                   !* Fill boundary data arrays with AGRIF 
     459      !                                   ! ------------------------------------- 
    469460#if defined key_agrif 
    470461         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) 

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

    r11239 r15764  
    2222   USE obs_read_surf            ! Reading and allocation of surface obs 
    2323   USE obs_readmdt              ! Reading and allocation of MDT for SLA. 
     24   USE obs_readsnowdepth        ! Get model snow depth for conversion of freeboard to ice thickness 
    2425   USE obs_prep                 ! Preparation of obs. (grid search etc). 
    2526   USE obs_oper                 ! Observation operators 
    … …  
    3334   USE mpp_map                  ! MPP mapping 
    3435   USE lib_mpp                  ! For ctl_warn/stop 
     36   USE tradmp                   ! For climatological temperature & salinity 
    3537 
    3638   IMPLICIT NONE 
    … …  
    4446 
    4547   !! * Module variables 
    46    LOGICAL, PUBLIC :: & 
    47       &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6. 
    48    LOGICAL :: ln_diaobs            !: Logical switch for the obs operator 
    49    LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
    50    LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
    51    LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres 
    52    LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
    53    LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
    54    LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
    55  
    56    REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
    57    REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 
    58    REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint 
    59    REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint 
    60    REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint 
    61    REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint 
    62    REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
    63    REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
    64    REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
    65    REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
    66  
    67    INTEGER :: nn_1dint         !: Vertical interpolation method 
    68    INTEGER :: nn_2dint_default !: Default horizontal interpolation method 
    69    INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default) 
    70    INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
    71    INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
    72    INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default) 
    73   
     48   LOGICAL, PUBLIC :: &  
     49      &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6.  
     50   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator  
     51   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs  
     52   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres  
     53   LOGICAL :: ln_sla_fp_indegs     !: T=> SLA obs footprint size specified in degrees, F=> in metres  
     54   LOGICAL :: ln_sst_fp_indegs     !: T=> SST obs footprint size specified in degrees, F=> in metres  
     55   LOGICAL :: ln_sss_fp_indegs     !: T=> SSS obs footprint size specified in degrees, F=> in metres  
     56   LOGICAL :: ln_ssv_fp_indegs     !: T=> SSV obs footprint size specified in degrees, F=> in metres     
     57   LOGICAL :: ln_sic_fp_indegs     !: T=> SIC obs footprint size specified in degrees, F=> in metres  
     58   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres  
     59   LOGICAL :: ln_fbd_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres  
     60   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology  
     61   LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal  
     62 
     63   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint  
     64   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint  
     65   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint  
     66   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint  
     67   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint  
     68   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint  
     69   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint  
     70   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint  
     71   REAL(wp) :: rn_ssv_avglamscl     !: E/W diameter of SSV observation footprint  
     72   REAL(wp) :: rn_ssv_avgphiscl     !: N/S diameter of SSV observation footprint  
     73   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of SIC observation footprint  
     74   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of SIC observation footprint  
     75   REAL(wp) :: rn_sit_avglamscl     !: E/W diameter of SIT observation footprint  
     76   REAL(wp) :: rn_sit_avgphiscl     !: N/S diameter of SIT observation footprint  
     77   REAL(wp) :: rn_fbd_avglamscl     !: E/W diameter of FBD observation footprint  
     78   REAL(wp) :: rn_fbd_avgphiscl     !: N/S diameter of FBD observation footprint  
     79   REAL(wp), PUBLIC :: &  
     80      &        MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data.  
     81 
     82 
     83   INTEGER :: nn_1dint         !: Vertical interpolation method  
     84   INTEGER :: nn_2dint_default !: Default horizontal interpolation method  
     85   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default)  
     86   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default)  
     87   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default)  
     88   INTEGER :: nn_2dint_ssv     !: SSV horizontal interpolation method (-1 = default)     
     89   INTEGER :: nn_2dint_sic     !: SIC horizontal interpolation method (-1 = default)  
     90   INTEGER :: nn_2dint_sit     !: SIT horizontal interpolation method (-1 = default)  
     91   INTEGER :: nn_2dint_fbd     !: FBD horizontal interpolation method (-1 = default) 
     92 
    7493   INTEGER, DIMENSION(imaxavtypes) :: & 
    7594      & nn_profdavtypes      !: Profile data types representing a daily average 
    … …  
    130149      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    131150      !!---------------------------------------------------------------------- 
     151#if defined key_cice  
     152USE sbc_oce, ONLY : &        ! CICE variables  
     153   & thick_s                 ! snow depth for freeboard conversion   
     154#endif 
    132155 
    133156      IMPLICIT NONE 
    134157 
    135158      !! * Local declarations 
    136       INTEGER, PARAMETER :: & 
    137          & jpmaxnfiles = 1000    ! Maximum number of files for each obs type 
    138       INTEGER, DIMENSION(:), ALLOCATABLE :: & 
    139          & ifilesprof, &         ! Number of profile files 
    140          & ifilessurf            ! Number of surface files 
    141       INTEGER :: ios             ! Local integer output status for namelist read 
    142       INTEGER :: jtype           ! Counter for obs types 
    143       INTEGER :: jvar            ! Counter for variables 
    144       INTEGER :: jfile           ! Counter for files 
    145       INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
    146       INTEGER :: n2dint_type     ! Local version of nn_2dint* 
    147  
    148       CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
    149          & cn_profbfiles,      & ! T/S profile input filenames 
    150          & cn_sstfbfiles,      & ! Sea surface temperature input filenames 
    151          & cn_slafbfiles,      & ! Sea level anomaly input filenames 
    152          & cn_sicfbfiles,      & ! Seaice concentration input filenames 
    153          & cn_velfbfiles,      & ! Velocity profile input filenames 
    154          & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
    155          & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames 
    156          & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames 
    157          & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames 
    158          & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames 
    159          & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 
    160          & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames 
    161          & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames 
    162          & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames 
    163          & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames 
    164          & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames 
    165          & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 
    166          & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames 
    167          & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames 
    168          & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames 
    169          & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 
    170          & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames 
    171          & cn_pno3fbfiles,     & ! Profile nitrate input filenames 
    172          & cn_psi4fbfiles,     & ! Profile silicate input filenames 
    173          & cn_ppo4fbfiles,     & ! Profile phosphate input filenames 
    174          & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames 
    175          & cn_palkfbfiles,     & ! Profile alkalinity input filenames 
    176          & cn_pphfbfiles,      & ! Profile pH input filenames 
    177          & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames 
    178          & cn_sstbiasfiles       ! SST bias input filenames 
    179  
    180       CHARACTER(LEN=128) :: & 
    181          & cn_altbiasfile        ! Altimeter bias input filename 
    182  
    183  
    184       LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
    185       LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
    186       LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
    187       LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
    188       LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
    189       LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
    190       LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
    191       LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs 
    192       LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs 
    193       LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs 
    194       LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs 
    195       LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 
    196       LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs 
    197       LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs 
    198       LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs 
    199       LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs 
    200       LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs 
    201       LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 
    202       LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs 
    203       LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs 
    204       LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs 
    205       LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs 
    206       LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs 
    207       LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs 
    208       LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs 
    209       LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs 
    210       LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs 
    211       LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs 
    212       LOGICAL :: ln_pph          ! Logical switch for profile pH obs 
    213       LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs 
    214       LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
    215       LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
    216       LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
    217       LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
    218       LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
    219       LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
    220  
    221       REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
    222       REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
    223  
    224       REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
    225       REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
    226  
    227       CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
    228          & clproffiles, &        ! Profile filenames 
    229          & clsurffiles           ! Surface filenames 
    230  
    231       LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read 
    232       LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 
    233       LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
    234  
    235       REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    236          & zglam                 ! Model longitudes for profile variables 
    237       REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    238          & zgphi                 ! Model latitudes for profile variables 
    239       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    240          & zmask                 ! Model land/sea mask associated with variables 
    241  
    242  
    243       NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    244          &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
    245          &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
    246          &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
    247          &            ln_slchlpic, ln_schltot,                        & 
    248          &            ln_slphytot, ln_slphydia, ln_slphynon,          & 
    249          &            ln_sspm,     ln_sfco2,    ln_spco2,             & 
    250          &            ln_plchltot, ln_pchltot,  ln_pno3,              & 
    251          &            ln_psi4,     ln_ppo4,     ln_pdic,              & 
    252          &            ln_palk,     ln_pph,      ln_po2,               & 
    253          &            ln_altbias, ln_sstbias, ln_nea,                 & 
    254          &            ln_grid_global, ln_grid_search_lookup,          & 
    255          &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    256          &            ln_sstnight, ln_default_fp_indegs,              & 
    257          &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    258          &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    259          &            cn_profbfiles, cn_slafbfiles,                   & 
    260          &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    261          &            cn_velfbfiles, cn_sssfbfiles,                   & 
    262          &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
    263          &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         & 
    264          &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         & 
    265          &            cn_slchlpicfbfiles, cn_schltotfbfiles,          & 
    266          &            cn_slphytotfbfiles, cn_slphydiafbfiles,         & 
    267          &            cn_slphynonfbfiles, cn_sspmfbfiles,             & 
    268          &            cn_sfco2fbfiles, cn_spco2fbfiles,               & 
    269          &            cn_plchltotfbfiles, cn_pchltotfbfiles,          & 
    270          &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 
    271          &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  & 
    272          &            cn_po2fbfiles,                                  & 
    273          &            cn_sstbiasfiles, cn_altbiasfile,                & 
    274          &            cn_gridsearchfile, rn_gridsearchres,            & 
    275          &            rn_dobsini, rn_dobsend,                         & 
    276          &            rn_default_avglamscl, rn_default_avgphiscl,     & 
    277          &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
    278          &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
    279          &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
    280          &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
    281          &            nn_1dint, nn_2dint_default,                     & 
    282          &            nn_2dint_sla, nn_2dint_sst,                     & 
    283          &            nn_2dint_sss, nn_2dint_sic,                     & 
    284          &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
     159      INTEGER, PARAMETER :: &  
     160         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type  
     161      INTEGER, DIMENSION(:), ALLOCATABLE :: &  
     162         & ifilesprof, &         ! Number of profile files  
     163         & ifilessurf            ! Number of surface files  
     164      INTEGER :: ios             ! Local integer output status for namelist read  
     165      INTEGER :: jtype           ! Counter for obs types  
     166      INTEGER :: jvar            ! Counter for variables  
     167      INTEGER :: jfile           ! Counter for files  
     168      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply  
     169      INTEGER :: n2dint_type     ! Local version of nn_2dint*  
     170 
     171      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &  
     172         & cn_profbfiles,      & ! T/S profile input filenames  
     173         & cn_sstfbfiles,      & ! Sea surface temperature input filenames  
     174         & cn_slafbfiles,      & ! Sea level anomaly input filenames  
     175         & cn_sicfbfiles,      & ! Seaice concentration input filenames  
     176         & cn_sitfbfiles,      & ! Seaice thickness input filenames  
     177         & cn_fbdfbfiles,      & ! Seaice freeboard input filenames  
     178         & cn_velfbfiles,      & ! Velocity profile input filenames  
     179         & cn_sssfbfiles,      & ! Sea surface salinity input filenames  
     180         & cn_ssvfbfiles,      & ! Sea surface velocity input filenames           
     181         & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames  
     182         & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames  
     183         & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames  
     184         & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames  
     185         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames  
     186         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames  
     187         & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames  
     188         & cn_schltotfbfiles,  & ! Surface total chlorophyll        input filenames  
     189         & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames  
     190         & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames  
     191         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames  
     192         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames  
     193         & cn_skd490fbfiles,   & ! Surface Kd490 input filenames  
     194         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames  
     195         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames  
     196         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames  
     197         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames  
     198         & cn_pno3fbfiles,     & ! Profile nitrate input filenames  
     199         & cn_psi4fbfiles,     & ! Profile silicate input filenames  
     200         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames  
     201         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames  
     202         & cn_palkfbfiles,     & ! Profile alkalinity input filenames  
     203         & cn_pphfbfiles,      & ! Profile pH input filenames  
     204         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames  
     205         & cn_sstbiasfiles       ! SST bias input filenames  
     206 
     207      CHARACTER(LEN=128) :: &  
     208         & cn_altbiasfile        ! Altimeter bias input filename  
     209 
     210 
     211      LOGICAL :: ln_seaicetypes = .FALSE.          ! Logical switch indicating data type is sea ice  
     212      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles  
     213      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles  
     214      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies   
     215      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature  
     216      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration  
     217      LOGICAL :: ln_sit          ! Logical switch for sea ice thickness  
     218      LOGICAL :: ln_fbd          ! Logical switch for sea ice freeboard  
     219      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs  
     220      LOGICAL :: ln_ssv          ! Logical switch for sea surface velocity obs        
     221      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs  
     222      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs  
     223      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs  
     224      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs  
     225      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs  
     226      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs  
     227      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs  
     228      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs  
     229      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs  
     230      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs  
     231      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs  
     232      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs  
     233      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs  
     234      LOGICAL :: ln_skd490       ! Logical switch for surface Kd490  
     235      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs  
     236      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs  
     237      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs  
     238      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs  
     239      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs  
     240      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs  
     241      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs  
     242      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs  
     243      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs  
     244      LOGICAL :: ln_pph          ! Logical switch for profile pH obs  
     245      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs  
     246      LOGICAL :: ln_nea          ! Logical switch to remove obs near land  
     247      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias  
     248      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST  
     249      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files  
     250      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs  
     251      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary  
     252 
     253      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS  
     254      REAL(dp) :: rn_dobsend     ! Obs window end date YYYYMMDD.HHMMSS  
     255 
     256      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl  
     257      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl  
     258 
     259      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &  
     260         & clproffiles, &        ! Profile filenames  
     261         & clsurffiles           ! Surface filenames  
     262      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars  ! Expected variable names  
     263 
     264      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read  
     265      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs  
     266      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)  
     267      LOGICAL :: ltype_clim      ! Local version of ln_output_clim  
     268 
     269      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     270         & zglam                 ! Model longitudes for profile variables  
     271      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     272         & zgphi                 ! Model latitudes for profile variables  
     273      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &  
     274         & zmask                 ! Model land/sea mask associated with variables  
     275      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     276         & zmask_surf            ! Surface model land/sea mask associated with variables  
     277 
     278 
     279      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &  
     280         &            ln_sst, ln_sic, ln_sit, ln_fbd,                 &  
     281         &            ln_sss, ln_ssv, ln_vel3d,                       &  
     282         &            ln_slchltot, ln_slchldia, ln_slchlnon,          &  
     283         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &  
     284         &            ln_slchlpic, ln_schltot,                        &  
     285         &            ln_slphytot, ln_slphydia, ln_slphynon,          &  
     286         &            ln_sspm,     ln_sfco2,    ln_spco2,             &  
     287         &            ln_skd490,                                      &  
     288         &            ln_plchltot, ln_pchltot,  ln_pno3,              &  
     289         &            ln_psi4,     ln_ppo4,     ln_pdic,              &  
     290         &            ln_palk,     ln_pph,      ln_po2,               &  
     291         &            ln_altbias, ln_sstbias, ln_nea,                 &  
     292         &            ln_grid_global, ln_grid_search_lookup,          &  
     293         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &  
     294         &            ln_sstnight,  ln_output_clim,                   &  
     295         &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     &  
     296         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &  
     297         &            ln_sss_fp_indegs, ln_ssv_fp_indegs,             &  
     298         &            ln_sic_fp_indegs, ln_sit_fp_indegs,             &  
     299         &            ln_fbd_fp_indegs,                               &  
     300         &            cn_profbfiles, cn_slafbfiles,                   &  
     301         &            cn_sstfbfiles, cn_sicfbfiles,                   &  
     302         &            cn_sitfbfiles, cn_fbdfbfiles,                   &  
     303         &            cn_velfbfiles, cn_sssfbfiles, cn_ssvfbfiles,    &  
     304         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &  
     305         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &  
     306         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         &  
     307         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          &  
     308         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         &  
     309         &            cn_slphynonfbfiles, cn_sspmfbfiles,             &  
     310         &            cn_skd490fbfiles,                               &  
     311         &            cn_sfco2fbfiles, cn_spco2fbfiles,               &  
     312         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          &  
     313         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, &  
     314         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  &  
     315         &            cn_po2fbfiles,                                  &  
     316         &            cn_sstbiasfiles, cn_altbiasfile,                &  
     317         &            cn_gridsearchfile, rn_gridsearchres,            &  
     318         &            rn_dobsini, rn_dobsend,                         &  
     319         &            rn_default_avglamscl, rn_default_avgphiscl,     &  
     320         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &  
     321         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &  
     322         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &  
     323         &            rn_ssv_avglamscl, rn_ssv_avgphiscl,             &           
     324         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &  
     325         &            rn_sit_avglamscl, rn_sit_avgphiscl,             &  
     326         &            rn_fbd_avglamscl, rn_fbd_avgphiscl,             &  
     327         &            nn_1dint, nn_2dint_default,                     &  
     328         &            nn_2dint_sla, nn_2dint_sst,                     &  
     329         &            nn_2dint_sss, nn_2dint_ssv,                     &  
     330         &            nn_2dint_sic, nn_2dint_sit,                     &  
     331         &            nn_2dint_fbd,                                   &  
     332         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &  
    285333         &            nn_profdavtypes 
    286334 
    … …  
    289337      !----------------------------------------------------------------------- 
    290338 
    291       ! Some namelist arrays need initialising 
    292       cn_profbfiles(:)      = '' 
    293       cn_slafbfiles(:)      = '' 
    294       cn_sstfbfiles(:)      = '' 
    295       cn_sicfbfiles(:)      = '' 
    296       cn_velfbfiles(:)      = '' 
    297       cn_sssfbfiles(:)      = '' 
    298       cn_slchltotfbfiles(:) = '' 
    299       cn_slchldiafbfiles(:) = '' 
    300       cn_slchlnonfbfiles(:) = '' 
    301       cn_slchldinfbfiles(:) = '' 
    302       cn_slchlmicfbfiles(:) = '' 
    303       cn_slchlnanfbfiles(:) = '' 
    304       cn_slchlpicfbfiles(:) = '' 
    305       cn_schltotfbfiles(:)  = '' 
    306       cn_slphytotfbfiles(:) = '' 
    307       cn_slphydiafbfiles(:) = '' 
    308       cn_slphynonfbfiles(:) = '' 
    309       cn_sspmfbfiles(:)     = '' 
    310       cn_sfco2fbfiles(:)    = '' 
    311       cn_spco2fbfiles(:)    = '' 
    312       cn_plchltotfbfiles(:) = '' 
    313       cn_pchltotfbfiles(:)  = '' 
    314       cn_pno3fbfiles(:)     = '' 
    315       cn_psi4fbfiles(:)     = '' 
    316       cn_ppo4fbfiles(:)     = '' 
    317       cn_pdicfbfiles(:)     = '' 
    318       cn_palkfbfiles(:)     = '' 
    319       cn_pphfbfiles(:)      = '' 
    320       cn_po2fbfiles(:)      = '' 
    321       cn_sstbiasfiles(:)    = '' 
    322       nn_profdavtypes(:)    = -1 
    323  
    324       CALL ini_date( rn_dobsini ) 
    325       CALL fin_date( rn_dobsend ) 
    326  
    327       ! Read namelist namobs : control observation diagnostics 
     339      cn_profbfiles(:)      = ''  
     340      cn_slafbfiles(:)      = ''  
     341      cn_sstfbfiles(:)      = ''  
     342      cn_sicfbfiles(:)      = ''  
     343      cn_sitfbfiles(:)      = ''  
     344      cn_fbdfbfiles(:)      = ''  
     345      cn_velfbfiles(:)      = ''  
     346      cn_sssfbfiles(:)      = ''  
     347      cn_ssvfbfiles(:)      = ''        
     348      cn_slchltotfbfiles(:) = ''  
     349      cn_slchldiafbfiles(:) = ''  
     350      cn_slchlnonfbfiles(:) = ''  
     351      cn_slchldinfbfiles(:) = ''  
     352      cn_slchlmicfbfiles(:) = ''  
     353      cn_slchlnanfbfiles(:) = ''  
     354      cn_slchlpicfbfiles(:) = ''  
     355      cn_schltotfbfiles(:)  = ''  
     356      cn_slphytotfbfiles(:) = ''  
     357      cn_slphydiafbfiles(:) = ''  
     358      cn_slphynonfbfiles(:) = ''  
     359      cn_sspmfbfiles(:)     = ''  
     360      cn_skd490fbfiles(:)   = ''  
     361      cn_sfco2fbfiles(:)    = ''  
     362      cn_spco2fbfiles(:)    = ''  
     363      cn_plchltotfbfiles(:) = ''  
     364      cn_pchltotfbfiles(:)  = ''  
     365      cn_pno3fbfiles(:)     = ''  
     366      cn_psi4fbfiles(:)     = ''  
     367      cn_ppo4fbfiles(:)     = ''  
     368      cn_pdicfbfiles(:)     = ''  
     369      cn_palkfbfiles(:)     = ''  
     370      cn_pphfbfiles(:)      = ''  
     371      cn_po2fbfiles(:)      = ''  
     372      cn_sstbiasfiles(:)    = ''  
     373      nn_profdavtypes(:)    = -1  
     374 
     375      CALL ini_date( rn_dobsini )  
     376      CALL fin_date( rn_dobsend )  
     377 
     378      ! Read namelist namobs : control observation diagnostics  
    328379      REWIND( numnam_ref )   ! Namelist namobs in reference namelist 
    329380      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 
    … …  
    351402         WRITE(numout,*) '~~~~~~~~~~~~' 
    352403         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters'  
    353          WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d 
    354          WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d 
    355          WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
    356          WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
    357          WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    358          WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    359          WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
    360          WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot 
    361          WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia 
    362          WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon 
    363          WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin 
    364          WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic 
    365          WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan 
    366          WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic 
    367          WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot 
    368          WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot 
    369          WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia 
    370          WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 
    371          WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm 
    372          WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2 
    373          WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2 
    374          WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot 
    375          WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot 
    376          WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3 
    377          WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4 
    378          WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4 
    379          WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic 
    380          WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk 
    381          WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph 
    382          WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2 
    383          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     404         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d  
     405         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d  
     406         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla  
     407         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst  
     408         WRITE(numout,*) '             Logical switch for SIC observations                      ln_sic = ', ln_sic  
     409         WRITE(numout,*) '             Logical switch for SIT observations                      ln_sit = ', ln_sit  
     410         WRITE(numout,*) '             Logical switch for FBD observations                      ln_fbd = ', ln_fbd  
     411         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d  
     412         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss  
     413         WRITE(numout,*) '             Logical switch for SSV observations                      ln_ssv = ', ln_ssv           
     414         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot  
     415         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia  
     416         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon  
     417         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin  
     418         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic  
     419         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan  
     420         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic  
     421         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot  
     422         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot  
     423         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia  
     424         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon  
     425         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm  
     426         WRITE(numout,*) '             Logical switch for surface Kd490 observations         ln_skd490 = ', ln_skd490  
     427         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2  
     428         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2  
     429         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot  
     430         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot  
     431         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3  
     432         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4  
     433         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4  
     434         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic  
     435         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk  
     436         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph  
     437         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2  
     438         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global  
    384439         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    385440         IF (ln_grid_search_lookup) & 
    386441            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
    387          WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini 
    388          WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
    389          WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
    390          WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
    391          WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
    392          WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
    393          WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
    394          WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
    395          WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
    396          WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
    397          WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs 
    398          WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl 
    399          WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl 
    400          WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs 
    401          WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl 
    402          WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl 
    403          WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs 
    404          WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl 
    405          WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
    406          WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
    407          WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
    408          WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    409          WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    410          WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    411          WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    412          WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
    413          WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    414          WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    415          WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    416          WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
     442            WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini  
     443            WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend  
     444            WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint  
     445            WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default  
     446            WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla  
     447            WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst  
     448            WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss  
     449            WRITE(numout,*) '             Type of horizontal interpolation method for SSV    nn_2dint_ssv = ', nn_2dint_ssv           
     450            WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic  
     451            WRITE(numout,*) '             Type of horizontal interpolation method for SIT    nn_2dint_sit = ', nn_2dint_sit  
     452            WRITE(numout,*) '             Type of horizontal interpolation method for FBD    nn_2dint_fbd = ', nn_2dint_fbd  
     453            WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl  
     454            WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl  
     455            WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs  
     456            WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl  
     457            WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl  
     458            WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs  
     459            WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl  
     460            WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl  
     461            WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs  
     462            WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl  
     463            WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl  
     464            WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs  
     465            WRITE(numout,*) '             SIT E/W diameter of obs footprint              rn_sit_avglamscl = ', rn_sit_avglamscl  
     466            WRITE(numout,*) '             SIT N/S diameter of obs footprint              rn_sit_avgphiscl = ', rn_sit_avgphiscl  
     467            WRITE(numout,*) '             SIT obs footprint in deg [T] or m [F]          ln_sit_fp_indegs = ', ln_sit_fp_indegs  
     468            WRITE(numout,*) '             FBD E/W diameter of obs footprint              rn_fbd_avglamscl = ', rn_fbd_avglamscl  
     469            WRITE(numout,*) '             FBD N/S diameter of obs footprint              rn_fbd_avgphiscl = ', rn_fbd_avgphiscl  
     470            WRITE(numout,*) '             FBD obs footprint in deg [T] or m [F]          ln_fbd_fp_indegs = ', ln_fbd_fp_indegs  
     471            WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea  
     472            WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject  
     473            WRITE(numout,*) '             MSSH correction scheme nn_msshc = ', nn_msshc  
     474            WRITE(numout,*) '             MDT  correction rn_mdtcorr = ', rn_mdtcorr  
     475            WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff  
     476            WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias  
     477            WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias  
     478            WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis  
     479            WRITE(numout,*) '             Daily average types nn_profdavtypes = ', nn_profdavtypes  
     480            WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight  
     481            WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim  
     482            WRITE(numout,*) '             Logical switch for time-mean of SLA        ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg 
    417483      ENDIF 
    418484      !----------------------------------------------------------------------- 
    … …  
    421487      !----------------------------------------------------------------------- 
    422488 
    423       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          & 
    424          &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
    425          &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
    426       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
    427          &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
    428          &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
    429          &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     & 
    430          &                  ln_sfco2,    ln_spco2 /) ) 
     489      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          &  
     490         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     &  
     491         &                  ln_pdic,     ln_palk,     ln_pph, ln_po2 /) )  
     492      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_fbd,             &  
     493         &                  ln_sss, ln_ssv,                                     &  
     494         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, &  
     495         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  &  
     496         &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     &  
     497         &                  ln_skd490,   ln_sfco2,    ln_spco2 /) )  
    431498 
    432499      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    … …  
    439506      ENDIF 
    440507 
     508      IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN  
     509         IF(lwp) WRITE(numout,cform_war)  
     510         IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', &  
     511            &                    ' so climatological T/S not available and will not be output'  
     512         nwarn = nwarn + 1  
     513         ln_output_clim = .FALSE.  
     514      ENDIF 
     515 
    441516      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    442517      IF ( nproftypes > 0 ) THEN 
    … …  
    535610            clsurffiles(jtype,:) = cn_sicfbfiles 
    536611         ENDIF 
     612         IF (ln_sit) THEN  
     613            jtype = jtype + 1  
     614            cobstypessurf(jtype) = 'sit'  
     615            clsurffiles(jtype,:) = cn_sitfbfiles  
     616         ENDIF  
     617         IF (ln_fbd) THEN  
     618            jtype = jtype + 1  
     619            cobstypessurf(jtype) = 'fbd'  
     620            clsurffiles(jtype,:) = cn_fbdfbfiles  
     621         ENDIF 
    537622         IF (ln_sss) THEN 
    538623            jtype = jtype + 1 
    … …  
    540625            clsurffiles(jtype,:) = cn_sssfbfiles 
    541626         ENDIF 
     627         IF (ln_ssv) THEN  
     628            jtype = jtype + 1  
     629            cobstypessurf(jtype) = 'ssv'  
     630            clsurffiles(jtype,:) = cn_ssvfbfiles  
     631         ENDIF 
    542632         IF (ln_slchltot) THEN 
    543633            jtype = jtype + 1 
    … …  
    599689            cobstypessurf(jtype) = 'sspm' 
    600690            clsurffiles(jtype,:) = cn_sspmfbfiles 
     691         ENDIF 
     692         IF (ln_skd490) THEN  
     693            jtype = jtype + 1  
     694            cobstypessurf(jtype) = 'skd490'  
     695            clsurffiles(jtype,:) = cn_skd490fbfiles  
    601696         ENDIF 
    602697         IF (ln_sfco2) THEN 
    … …  
    645740               ltype_fp_indegs = ln_sic_fp_indegs 
    646741               ltype_night     = .FALSE. 
     742            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN  
     743               IF ( nn_2dint_sit == -1 ) THEN  
     744                  n2dint_type  = nn_2dint_default  
     745               ELSE  
     746                  n2dint_type  = nn_2dint_sit  
     747               ENDIF  
     748               ztype_avglamscl = rn_sit_avglamscl  
     749               ztype_avgphiscl = rn_sit_avgphiscl  
     750               ltype_fp_indegs = ln_sit_fp_indegs  
     751               ltype_night     = .FALSE.  
     752            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN  
     753               IF ( nn_2dint_fbd == -1 ) THEN  
     754                  n2dint_type  = nn_2dint_default  
     755               ELSE  
     756                  n2dint_type  = nn_2dint_fbd  
     757               ENDIF  
     758               ztype_avglamscl = rn_fbd_avglamscl  
     759               ztype_avgphiscl = rn_fbd_avgphiscl  
     760               ltype_fp_indegs = ln_fbd_fp_indegs  
     761               ltype_night     = .FALSE. 
    647762            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
    648763               IF ( nn_2dint_sss == -1 ) THEN 
    … …  
    654769               ztype_avgphiscl = rn_sss_avgphiscl 
    655770               ltype_fp_indegs = ln_sss_fp_indegs 
     771               ltype_night     = .FALSE. 
     772            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN  
     773               IF ( nn_2dint_ssv == -1 ) THEN  
     774                  n2dint_type  = nn_2dint_default  
     775               ELSE  
     776                  n2dint_type  = nn_2dint_ssv  
     777               ENDIF  
     778               ztype_avglamscl = rn_ssv_avglamscl  
     779               ztype_avgphiscl = rn_ssv_avgphiscl  
     780               ltype_fp_indegs = ln_ssv_fp_indegs  
    656781               ltype_night     = .FALSE. 
    657782            ELSE 
    … …  
    715840         DO jtype = 1, nproftypes 
    716841 
     842            ltype_clim = .FALSE. 
     843 
    717844            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
    718845               nvarsprof(jtype) = 2 
    719846               nextrprof(jtype) = 1 
     847               IF ( ln_output_clim ) ltype_clim = .TRUE. 
    720848               ALLOCATE(llvar(nvarsprof(jtype))) 
     849               ALLOCATE(clvars(nvarsprof(jtype))) 
    721850               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
    722851               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
    … …  
    734863               nextrprof(jtype) = 2 
    735864               ALLOCATE(llvar(nvarsprof(jtype))) 
     865               ALLOCATE(clvars(nvarsprof(jtype))) 
    736866               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
    737867               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
    … …  
    739869               llvar(1)       = ln_vel3d 
    740870               llvar(2)       = ln_vel3d 
     871               clvars(1)      = 'UVEL'  
     872               clvars(2)      = 'VVEL' 
    741873               zglam(:,:,1)   = glamu(:,:) 
    742874               zglam(:,:,2)   = glamv(:,:) 
    … …  
    749881               nextrprof(jtype) = 0 
    750882               ALLOCATE(llvar(nvarsprof(jtype))) 
     883               ALLOCATE(clvars(nvarsprof(jtype))) 
    751884               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
    752885               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
    … …  
    756889               zgphi(:,:,1)   = gphit(:,:) 
    757890               zmask(:,:,:,1) = tmask(:,:,:) 
     891               IF ( TRIM(cobstypesprof(jtype)) == 'plchltot' ) THEN  
     892                  clvars(1) = 'PLCHLTOT'  
     893               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pchltot' )  THEN  
     894                  clvars(1) = 'PCHLTOT'  
     895               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pno3' ) THEN  
     896                  clvars(1) = 'PNO3'  
     897               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'psi4' ) THEN  
     898                  clvars(1) = 'PSI4'  
     899               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'ppo4' ) THEN  
     900                  clvars(1) = 'PPO4'  
     901               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pdic' ) THEN  
     902                  clvars(1) = 'PDIC'  
     903               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'palk' ) THEN  
     904                  clvars(1) = 'PALK'  
     905               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pph' ) THEN  
     906                  clvars(1) = 'PPH'  
     907               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'po2' ) THEN  
     908                  clvars(1) = 'PO2'  
     909               ENDIF 
    758910            ENDIF 
    759911 
    … …  
    763915               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
    764916               &               rn_dobsini, rn_dobsend, llvar, & 
    765                &               ln_ignmis, ln_s_at_t, .FALSE., & 
     917               &               ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, clvars, & 
    766918               &               kdailyavtypes = nn_profdavtypes ) 
    767919 
    … …  
    797949         DO jtype = 1, nsurftypes 
    798950 
    799             nvarssurf(jtype) = 1 
    800             nextrsurf(jtype) = 0 
    801             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
     951            ltype_clim = .FALSE.  
     952            IF ( ln_output_clim .AND. &  
     953               & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. &  
     954               &   ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) &  
     955               & ltype_clim = .TRUE.  
     956 
     957            IF ( (TRIM(cobstypessurf(jtype)) == 'sla') .OR. &  
     958               & (TRIM(cobstypessurf(jtype)) == 'sit') .OR. &  
     959               & (TRIM(cobstypessurf(jtype)) == 'fbd') ) THEN  
     960               nvarssurf(jtype) = 1  
     961               nextrsurf(jtype) = 2  
     962            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN  
     963               nvarssurf(jtype) = 2  
     964               nextrsurf(jtype) = 0  
     965            ELSE  
     966               nvarssurf(jtype) = 1  
     967               nextrsurf(jtype) = 0  
     968            ENDIF  
     969 
     970            ALLOCATE( clvars( nvarssurf(jtype) ) )  
     971            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zglam )  
     972            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zgphi )  
     973            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zmask_surf )  
     974   
     975            IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN  
     976               zglam(:,:,1) = glamu(:,:)  
     977               zglam(:,:,2) = glamv(:,:)  
     978               zgphi(:,:,1) = gphiu(:,:)  
     979               zgphi(:,:,2) = gphiv(:,:)  
     980               zmask_surf(:,:,1) = umask(:,:,1)  
     981               zmask_surf(:,:,2) = vmask(:,:,1)  
     982            ELSE              
     983               DO jvar = 1, nvarssurf(jtype)  
     984                  zglam(:,:,jvar) = glamt(:,:)  
     985                  zgphi(:,:,jvar) = gphit(:,:)  
     986                  zmask_surf(:,:,jvar) = tmask(:,:,1)                                         
     987               END DO 
     988            ENDIF 
     989 
     990            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN  
     991               clvars(1) = 'SLA'  
     992            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN  
     993               clvars(1) = 'SST'  
     994            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN  
     995               clvars(1) = 'ICECONC'  
     996            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN  
     997               clvars(1) = 'FBD'  
     998               ln_seaicetypes = .TRUE.  
     999            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN  
     1000               clvars(1) = 'SIT'  
     1001               ln_seaicetypes = .TRUE.  
     1002            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN  
     1003               clvars(1) = 'SSS'  
     1004            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN  
     1005               clvars(1) = 'UVEL'  
     1006               clvars(2) = 'VVEL'             
     1007            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN  
     1008               clvars(1) = 'SLCHLTOT'  
     1009            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldia' ) THEN  
     1010               clvars(1) = 'SLCHLDIA'  
     1011            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnon' ) THEN  
     1012               clvars(1) = 'SLCHLNON'  
     1013            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldin' ) THEN  
     1014               clvars(1) = 'SLCHLDIN'  
     1015            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlmic' ) THEN  
     1016               clvars(1) = 'SLCHLMIC'  
     1017            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnan' ) THEN  
     1018               clvars(1) = 'SLCHLNAN'  
     1019            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlpic' ) THEN  
     1020               clvars(1) = 'SLCHLPIC'  
     1021            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'schltot' ) THEN  
     1022               clvars(1) = 'SCHLTOT'  
     1023            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphytot' ) THEN  
     1024               clvars(1) = 'SLPHYTOT'  
     1025            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphydia' ) THEN  
     1026               clvars(1) = 'SLPHYDIA'  
     1027            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphynon' ) THEN  
     1028               clvars(1) = 'SLPHYNON'  
     1029            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sspm' ) THEN  
     1030               clvars(1) = 'SSPM'  
     1031            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'skd490' ) THEN  
     1032               clvars(1) = 'SKD490'  
     1033            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sfco2' ) THEN  
     1034               clvars(1) = 'SFCO2'  
     1035            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'spco2' ) THEN  
     1036               clvars(1) = 'SPCO2'  
     1037            ENDIF 
    8021038 
    8031039            !Read in surface obs types 
    804             CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
    805                &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    806                &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    807                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
    808  
    809             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    810  
    811             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    812                CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
    813                IF ( ln_altbias ) & 
    814                   & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    815             ENDIF 
    816  
    817             IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
    818                jnumsstbias = 0 
    819                DO jfile = 1, jpmaxnfiles 
    820                   IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
    821                      &  jnumsstbias = jnumsstbias + 1 
    822                END DO 
    823                IF ( jnumsstbias == 0 ) THEN 
    824                   CALL ctl_stop("ln_sstbias set but no bias files to read in")     
     1040            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &  
     1041               &               clsurffiles(jtype,1:ifilessurf(jtype)), &  
     1042               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &  
     1043               &               rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., &  
     1044               &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars )  
     1045  
     1046            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), &  
     1047               &               jpi, jpj, &              
     1048               &               zmask_surf, zglam, zgphi, &  
     1049               &               ln_nea, ln_bound_reject, ln_seaicetypes )  
     1050  
     1051            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN  
     1052               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )  
     1053               IF ( ln_altbias ) &  
     1054                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )  
     1055            ENDIF  
     1056  
     1057#if defined key_cice  
     1058            IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN  
     1059               CALL obs_rea_snowdepth( surfdataqc(jtype), n2dintsurf(jtype), thick_s(:,:) )  
     1060            ENDIF   
     1061#endif  
     1062  
     1063            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN  
     1064               jnumsstbias = 0  
     1065               DO jfile = 1, jpmaxnfiles  
     1066                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &  
     1067                     &  jnumsstbias = jnumsstbias + 1  
     1068               END DO  
     1069               IF ( jnumsstbias == 0 ) THEN  
     1070                  CALL ctl_stop("ln_sstbias set but no bias files to read in")      
     1071  
    8251072               ENDIF 
    8261073 
    827                CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
    828                   &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
    829  
    830             ENDIF 
    831  
    832          END DO 
    833  
     1074               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &   
     1075                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )   
     1076  
     1077            ENDIF  
     1078              
     1079            DEALLOCATE( clvars )  
     1080            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zglam )  
     1081            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zgphi )  
     1082            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zmask_surf )  
     1083  
     1084         END DO  
     1085  
    8341086         DEALLOCATE( ifilessurf, clsurffiles ) 
    8351087 
    … …  
    8791131         & frld 
    8801132#endif 
    881 #if defined key_cice 
    882       USE sbc_oce, ONLY : fr_i     ! ice fraction 
    883 #endif 
    884 #if defined key_top 
    885       USE trc, ONLY :  &           ! Biogeochemical state variables 
    886          & trn 
    887 #endif 
    888 #if defined key_hadocc 
    889       USE par_hadocc               ! HadOCC parameters 
    890       USE trc, ONLY :  & 
    891          & HADOCC_CHL, & 
    892          & HADOCC_FCO2, & 
    893          & HADOCC_PCO2, & 
    894          & HADOCC_FILL_FLT 
    895       USE had_bgc_const, ONLY: c2n_p 
    896 #elif defined key_medusa 
    897       USE par_medusa               ! MEDUSA parameters 
    898       USE sms_medusa, ONLY: & 
    899          & xthetapn, & 
    900          & xthetapd 
    901 #if defined key_roam 
    902       USE sms_medusa, ONLY: & 
    903          & f2_pco2w, & 
    904          & f2_fco2w, & 
    905          & f3_pH 
    906 #endif 
    907 #elif defined key_fabm 
    908       USE par_fabm                 ! FABM parameters 
    909       USE fabm, ONLY: & 
    910          & fabm_get_interior_diagnostic_data 
    911 #endif 
    912 #if defined key_spm 
    913       USE par_spm, ONLY: &         ! Sediment parameters 
    914          & jp_spm 
     1133#if defined key_cice  
     1134      USE sbc_oce, ONLY : &        ! CICE variables  
     1135         & fr_i,          &        ! ice fraction  
     1136         & thick_i                 ! ice thickness  
     1137#endif  
     1138#if defined key_top  
     1139      USE trc, ONLY :  &           ! Biogeochemical state variables  
     1140         & trn  
     1141#endif  
     1142#if defined key_hadocc  
     1143      USE par_hadocc               ! HadOCC parameters  
     1144      USE trc, ONLY :  &  
     1145         & HADOCC_CHL, &  
     1146         & HADOCC_FCO2, &  
     1147         & HADOCC_PCO2, &  
     1148         & HADOCC_FILL_FLT  
     1149      USE had_bgc_const, ONLY: c2n_p  
     1150#elif defined key_medusa  
     1151      USE par_medusa               ! MEDUSA parameters  
     1152      USE sms_medusa, ONLY: &  
     1153         & xthetapn, &  
     1154         & xthetapd  
     1155#if defined key_roam  
     1156      USE sms_medusa, ONLY: &  
     1157         & f2_pco2w, &  
     1158         & f2_fco2w, &  
     1159         & f3_pH  
     1160#endif  
     1161#elif defined key_fabm  
     1162      USE par_fabm                 ! FABM parameters  
     1163#endif  
     1164#if defined key_spm  
     1165      USE par_spm, ONLY: &         ! Sediment parameters  
     1166         & jp_spm  
    9151167#endif 
    9161168 
    … …  
    9201172      INTEGER, INTENT(IN) :: kstp  ! Current timestep 
    9211173      !! * Local declarations 
    922       INTEGER :: idaystp           ! Number of timesteps per day 
    923       INTEGER :: jtype             ! Data loop variable 
    924       INTEGER :: jvar              ! Variable number 
    925       INTEGER :: ji, jj, jk        ! Loop counters 
    926       REAL(wp) :: tiny             ! small number 
    927       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    928          & zprofvar                ! Model values for variables in a prof ob 
    929       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    930          & zprofmask               ! Mask associated with zprofvar 
    931       REAL(wp), POINTER, DIMENSION(:,:) :: & 
    932          & zsurfvar, &             ! Model values equivalent to surface ob. 
    933          & zsurfmask               ! Mask associated with surface variable 
    934       REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    935          & zglam,    &             ! Model longitudes for prof variables 
    936          & zgphi                   ! Model latitudes for prof variables 
    937       LOGICAL :: llog10            ! Perform log10 transform of variable 
    938 #if defined key_fabm 
    939       REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    940          & pco2_3d                 ! 3D pCO2 from FABM 
     1174      INTEGER :: idaystp           ! Number of timesteps per day  
     1175      INTEGER :: imeanstp          ! Number of timesteps for sla averaging  
     1176      INTEGER :: jtype             ! Data loop variable  
     1177      INTEGER :: jvar              ! Variable number  
     1178      INTEGER :: ji, jj, jk        ! Loop counters  
     1179      REAL(wp) :: tiny             ! small number  
     1180      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &  
     1181         & zprofvar, &             ! Model values for variables in a prof ob  
     1182         & zprofclim               ! Climatology values for variables in a prof ob  
     1183      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &  
     1184         & zprofmask               ! Mask associated with zprofvar  
     1185      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     1186         & zsurfvar, &             ! Model values equivalent to surface ob.  
     1187         & zsurfclim, &            ! Climatology values for variables in a surface ob.  
     1188         & zsurfmask               ! Mask associated with surface variable  
     1189      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     1190         & zglam,    &             ! Model longitudes  
     1191         & zgphi                   ! Model latitudes  
     1192      LOGICAL :: llog10            ! Perform log10 transform of variable  
     1193#if defined key_fabm  
     1194      REAL(wp), POINTER, DIMENSION(:,:,:) :: &  
     1195         & fabm_3d                 ! 3D variable from FABM  
    9411196#endif 
    9421197 
    … …  
    9541209      !----------------------------------------------------------------------- 
    9551210 
    956       IF ( nproftypes > 0 ) THEN 
    957  
    958          DO jtype = 1, nproftypes 
    959  
    960             ! Allocate local work arrays 
    961             CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
    962             CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
    963             CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
    964             CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
    965              
    966             ! Defaults which might change 
    967             DO jvar = 1, profdataqc(jtype)%nvar 
    968                zprofmask(:,:,:,jvar) = tmask(:,:,:) 
    969                zglam(:,:,jvar)       = glamt(:,:) 
    970                zgphi(:,:,jvar)       = gphit(:,:) 
    971             END DO 
    972  
    973             SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
    974  
    975             CASE('prof') 
    976                zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 
    977                zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 
    978  
    979             CASE('vel') 
    980                zprofvar(:,:,:,1) = un(:,:,:) 
    981                zprofvar(:,:,:,2) = vn(:,:,:) 
    982                zprofmask(:,:,:,1) = umask(:,:,:) 
    983                zprofmask(:,:,:,2) = vmask(:,:,:) 
    984                zglam(:,:,1) = glamu(:,:) 
    985                zglam(:,:,2) = glamv(:,:) 
    986                zgphi(:,:,1) = gphiu(:,:) 
    987                zgphi(:,:,2) = gphiv(:,:) 
    988  
    989             CASE('plchltot') 
    990 #if defined key_hadocc 
    991                ! Chlorophyll from HadOCC 
    992                zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
    993 #elif defined key_medusa 
    994                ! Add non-diatom and diatom chlorophyll from MEDUSA 
    995                zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
    996 #elif defined key_fabm 
    997                ! Add all chlorophyll groups from ERSEM 
    998                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
    999                   &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
    1000 #else 
    1001                CALL ctl_stop( ' Trying to run plchltot observation operator', & 
    1002                   &           ' but no biogeochemical model appears to have been defined' ) 
    1003 #endif 
    1004                ! Take the log10 where we can, otherwise exclude 
    1005                tiny = 1.0e-20 
    1006                WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 
    1007                   zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:)) 
    1008                ELSEWHERE 
    1009                   zprofvar(:,:,:,:)  = obfillflt 
    1010                   zprofmask(:,:,:,:) = 0 
    1011                END WHERE 
    1012                ! Mask out model below any excluded values, 
    1013                ! to avoid interpolation issues 
    1014                DO jvar = 1, profdataqc(jtype)%nvar 
    1015                  DO jj = 1, jpj 
    1016                     DO ji = 1, jpi 
    1017                        depth_loop: DO jk = 1, jpk 
    1018                           IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 
    1019                              zprofmask(ji,jj,jk:jpk,jvar) = 0 
    1020                              EXIT depth_loop 
    1021                           ENDIF 
    1022                        END DO depth_loop 
    1023                     END DO 
    1024                  END DO 
    1025               END DO 
    1026  
    1027             CASE('pchltot') 
    1028 #if defined key_hadocc 
    1029                ! Chlorophyll from HadOCC 
    1030                zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
    1031 #elif defined key_medusa 
    1032                ! Add non-diatom and diatom chlorophyll from MEDUSA 
    1033                zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
    1034 #elif defined key_fabm 
    1035                ! Add all chlorophyll groups from ERSEM 
    1036                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
    1037                   &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
    1038 #else 
    1039                CALL ctl_stop( ' Trying to run pchltot observation operator', & 
    1040                   &           ' but no biogeochemical model appears to have been defined' ) 
    1041 #endif 
    1042  
    1043             CASE('pno3') 
    1044 #if defined key_hadocc 
    1045                ! Dissolved inorganic nitrogen from HadOCC 
    1046                zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 
    1047 #elif defined key_medusa 
    1048                ! Dissolved inorganic nitrogen from MEDUSA 
    1049                zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 
    1050 #elif defined key_fabm 
    1051                ! Nitrate from ERSEM 
    1052                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 
    1053 #else 
    1054                CALL ctl_stop( ' Trying to run pno3 observation operator', & 
    1055                   &           ' but no biogeochemical model appears to have been defined' ) 
    1056 #endif 
    1057  
    1058             CASE('psi4') 
    1059 #if defined key_hadocc 
    1060                CALL ctl_stop( ' Trying to run psi4 observation operator', & 
    1061                   &           ' but HadOCC does not simulate silicate' ) 
    1062 #elif defined key_medusa 
    1063                ! Silicate from MEDUSA 
    1064                zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 
    1065 #elif defined key_fabm 
    1066                ! Silicate from ERSEM 
    1067                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 
    1068 #else 
    1069                CALL ctl_stop( ' Trying to run psi4 observation operator', & 
    1070                   &           ' but no biogeochemical model appears to have been defined' ) 
    1071 #endif 
    1072  
    1073             CASE('ppo4') 
    1074 #if defined key_hadocc 
    1075                CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
    1076                   &           ' but HadOCC does not simulate phosphate' ) 
    1077 #elif defined key_medusa 
    1078                CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
    1079                   &           ' but MEDUSA does not simulate phosphate' ) 
    1080 #elif defined key_fabm 
    1081                ! Phosphate from ERSEM 
    1082                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 
    1083 #else 
    1084                CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
    1085                   &           ' but no biogeochemical model appears to have been defined' ) 
    1086 #endif 
    1087  
    1088             CASE('pdic') 
    1089 #if defined key_hadocc 
    1090                ! Dissolved inorganic carbon from HadOCC 
    1091                zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 
    1092 #elif defined key_medusa 
    1093                ! Dissolved inorganic carbon from MEDUSA 
    1094                zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 
    1095 #elif defined key_fabm 
    1096                ! Dissolved inorganic carbon from ERSEM 
    1097                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 
    1098 #else 
    1099                CALL ctl_stop( ' Trying to run pdic observation operator', & 
    1100                   &           ' but no biogeochemical model appears to have been defined' ) 
    1101 #endif 
    1102  
    1103             CASE('palk') 
    1104 #if defined key_hadocc 
    1105                ! Alkalinity from HadOCC 
    1106                zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 
    1107 #elif defined key_medusa 
    1108                ! Alkalinity from MEDUSA 
    1109                zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 
    1110 #elif defined key_fabm 
    1111                ! Alkalinity from ERSEM 
    1112                zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 
    1113 #else 
    1114                CALL ctl_stop( ' Trying to run palk observation operator', & 
    1115                   &           ' but no biogeochemical model appears to have been defined' ) 
    1116 #endif 
    1117  
    1118             CASE('pph') 
    1119 #if defined key_hadocc 
    1120                CALL ctl_stop( ' Trying to run pph observation operator', & 
    1121                   &           ' but HadOCC has no pH diagnostic defined' ) 
    1122 #elif defined key_medusa && defined key_roam 
    1123                ! pH from MEDUSA 
    1124                zprofvar(:,:,:,1) = f3_pH(:,:,:) 
    1125 #elif defined key_fabm 
    1126                ! pH from ERSEM 
    1127                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 
    1128 #else 
    1129                CALL ctl_stop( ' Trying to run pph observation operator', & 
    1130                   &           ' but no biogeochemical model appears to have been defined' ) 
    1131 #endif 
    1132  
    1133             CASE('po2') 
    1134 #if defined key_hadocc 
    1135                CALL ctl_stop( ' Trying to run po2 observation operator', & 
    1136                   &           ' but HadOCC does not simulate oxygen' ) 
    1137 #elif defined key_medusa 
    1138                ! Oxygen from MEDUSA 
    1139                zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 
    1140 #elif defined key_fabm 
    1141                ! Oxygen from ERSEM 
    1142                zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 
    1143 #else 
    1144                CALL ctl_stop( ' Trying to run po2 observation operator', & 
    1145                   &           ' but no biogeochemical model appears to have been defined' ) 
    1146 #endif 
    1147  
    1148             CASE DEFAULT 
    1149                CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    1150  
    1151             END SELECT 
    1152  
    1153             DO jvar = 1, profdataqc(jtype)%nvar 
    1154                CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    1155                   &               nit000, idaystp, jvar,                   & 
    1156                   &               zprofvar(:,:,:,jvar),                    & 
    1157                   &               fsdept(:,:,:), fsdepw(:,:,:),            &  
    1158                   &               zprofmask(:,:,:,jvar),                   & 
    1159                   &               zglam(:,:,jvar), zgphi(:,:,jvar),        & 
    1160                   &               nn_1dint, nn_2dint_default,              & 
    1161                   &               kdailyavtypes = nn_profdavtypes ) 
    1162             END DO 
    1163  
    1164             CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
    1165             CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
    1166             CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
    1167             CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
    1168  
    1169          END DO 
    1170  
    1171       ENDIF 
    1172  
    1173       IF ( nsurftypes > 0 ) THEN 
    1174  
    1175          !Allocate local work arrays 
    1176          CALL wrk_alloc( jpi, jpj, zsurfvar ) 
    1177          CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    1178 #if defined key_fabm 
    1179          CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 
    1180 #endif 
    1181  
    1182          DO jtype = 1, nsurftypes 
    1183  
    1184             !Defaults which might be changed 
    1185             zsurfmask(:,:) = tmask(:,:,1) 
    1186             llog10 = .FALSE. 
    1187  
    1188             SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    1189             CASE('sst') 
    1190                zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    1191             CASE('sla') 
    1192                zsurfvar(:,:) = sshn(:,:) 
    1193             CASE('sss') 
    1194                zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    1195             CASE('sic') 
    1196                IF ( kstp == 0 ) THEN 
    1197                   IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
    1198                      CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
    1199                         &           'time-step but some obs are valid then.' ) 
    1200                      WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
    1201                         &           ' sea-ice obs will be missed' 
    1202                   ENDIF 
    1203                   surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
    1204                      &                        surfdataqc(jtype)%nsstp(1) 
    1205                   CYCLE 
    1206                ELSE 
    1207 #if defined key_cice 
    1208                   zsurfvar(:,:) = fr_i(:,:) 
    1209 #elif defined key_lim2 || defined key_lim3 
    1210                   zsurfvar(:,:) = 1._wp - frld(:,:) 
    1211 #else 
    1212                CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
    1213                   &           ' but no sea-ice model appears to have been defined' ) 
    1214 #endif 
    1215                ENDIF 
    1216  
    1217             CASE('slchltot') 
    1218 #if defined key_hadocc 
    1219                ! Surface chlorophyll from HadOCC 
    1220                zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
    1221 #elif defined key_medusa 
    1222                ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
    1223                zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
    1224 #elif defined key_fabm 
    1225                ! Add all surface chlorophyll groups from ERSEM 
    1226                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    1227                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    1228 #else 
    1229                CALL ctl_stop( ' Trying to run slchltot observation operator', & 
    1230                   &           ' but no biogeochemical model appears to have been defined' ) 
    1231 #endif 
    1232                llog10 = .TRUE. 
    1233  
    1234             CASE('slchldia') 
    1235 #if defined key_hadocc 
    1236                CALL ctl_stop( ' Trying to run slchldia observation operator', & 
    1237                   &           ' but HadOCC does not explicitly simulate diatoms' ) 
    1238 #elif defined key_medusa 
    1239                ! Diatom surface chlorophyll from MEDUSA 
    1240                zsurfvar(:,:) = trn(:,:,1,jpchd) 
    1241 #elif defined key_fabm 
    1242                ! Diatom surface chlorophyll from ERSEM 
    1243                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
    1244 #else 
    1245                CALL ctl_stop( ' Trying to run slchldia observation operator', & 
    1246                   &           ' but no biogeochemical model appears to have been defined' ) 
    1247 #endif 
    1248                llog10 = .TRUE. 
    1249  
    1250             CASE('slchlnon') 
    1251 #if defined key_hadocc 
    1252                CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
    1253                   &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
    1254 #elif defined key_medusa 
    1255                ! Non-diatom surface chlorophyll from MEDUSA 
    1256                zsurfvar(:,:) = trn(:,:,1,jpchn) 
    1257 #elif defined key_fabm 
    1258                ! Add all non-diatom surface chlorophyll groups from ERSEM 
    1259                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    1260                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    1261 #else 
    1262                CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
    1263                   &           ' but no biogeochemical model appears to have been defined' ) 
    1264 #endif 
    1265                llog10 = .TRUE. 
    1266  
    1267             CASE('slchldin') 
    1268 #if defined key_hadocc 
    1269                CALL ctl_stop( ' Trying to run slchldin observation operator', & 
    1270                   &           ' but HadOCC does not explicitly simulate dinoflagellates' ) 
    1271 #elif defined key_medusa 
    1272                CALL ctl_stop( ' Trying to run slchldin observation operator', & 
    1273                   &           ' but MEDUSA does not explicitly simulate dinoflagellates' ) 
    1274 #elif defined key_fabm 
    1275                ! Dinoflagellate surface chlorophyll from ERSEM 
    1276                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    1277 #else 
    1278                CALL ctl_stop( ' Trying to run slchldin observation operator', & 
    1279                   &           ' but no biogeochemical model appears to have been defined' ) 
    1280 #endif 
    1281                llog10 = .TRUE. 
    1282  
    1283             CASE('slchlmic') 
    1284 #if defined key_hadocc 
    1285                CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
    1286                   &           ' but HadOCC does not explicitly simulate microphytoplankton' ) 
    1287 #elif defined key_medusa 
    1288                CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
    1289                   &           ' but MEDUSA does not explicitly simulate microphytoplankton' ) 
    1290 #elif defined key_fabm 
    1291                ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 
    1292                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    1293 #else 
    1294                CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
    1295                   &           ' but no biogeochemical model appears to have been defined' ) 
    1296 #endif 
    1297                llog10 = .TRUE. 
    1298  
    1299             CASE('slchlnan') 
    1300 #if defined key_hadocc 
    1301                CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
    1302                   &           ' but HadOCC does not explicitly simulate nanophytoplankton' ) 
    1303 #elif defined key_medusa 
    1304                CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
    1305                   &           ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 
    1306 #elif defined key_fabm 
    1307                ! Nanophytoplankton surface chlorophyll from ERSEM 
    1308                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
    1309 #else 
    1310                CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
    1311                   &           ' but no biogeochemical model appears to have been defined' ) 
    1312 #endif 
    1313                llog10 = .TRUE. 
    1314  
    1315             CASE('slchlpic') 
    1316 #if defined key_hadocc 
    1317                CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
    1318                   &           ' but HadOCC does not explicitly simulate picophytoplankton' ) 
    1319 #elif defined key_medusa 
    1320                CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
    1321                   &           ' but MEDUSA does not explicitly simulate picophytoplankton' ) 
    1322 #elif defined key_fabm 
    1323                ! Picophytoplankton surface chlorophyll from ERSEM 
    1324                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
    1325 #else 
    1326                CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
    1327                   &           ' but no biogeochemical model appears to have been defined' ) 
    1328 #endif 
    1329                llog10 = .TRUE. 
    1330  
    1331             CASE('schltot') 
    1332 #if defined key_hadocc 
    1333                ! Surface chlorophyll from HadOCC 
    1334                zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
    1335 #elif defined key_medusa 
    1336                ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
    1337                zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
    1338 #elif defined key_fabm 
    1339                ! Add all surface chlorophyll groups from ERSEM 
    1340                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    1341                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    1342 #else 
    1343                CALL ctl_stop( ' Trying to run schltot observation operator', & 
    1344                   &           ' but no biogeochemical model appears to have been defined' ) 
    1345 #endif 
    1346  
    1347             CASE('slphytot') 
    1348 #if defined key_hadocc 
    1349                ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 
    1350                zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 
    1351 #elif defined key_medusa 
    1352                ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 
    1353                ! multiplied by C:N ratio for each 
    1354                zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
    1355 #elif defined key_fabm 
    1356                ! Add all surface phytoplankton carbon groups from ERSEM 
    1357                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
    1358                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
    1359 #else 
    1360                CALL ctl_stop( ' Trying to run slphytot observation operator', & 
    1361                   &           ' but no biogeochemical model appears to have been defined' ) 
    1362 #endif 
    1363                llog10 = .TRUE. 
    1364  
    1365             CASE('slphydia') 
    1366 #if defined key_hadocc 
    1367                CALL ctl_stop( ' Trying to run slphydia observation operator', & 
    1368                   &           ' but HadOCC does not explicitly simulate diatoms' ) 
    1369 #elif defined key_medusa 
    1370                ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
    1371                zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 
    1372 #elif defined key_fabm 
    1373                ! Diatom surface phytoplankton carbon from ERSEM 
    1374                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
    1375 #else 
    1376                CALL ctl_stop( ' Trying to run slphydia observation operator', & 
    1377                   &           ' but no biogeochemical model appears to have been defined' ) 
    1378 #endif 
    1379                llog10 = .TRUE. 
    1380  
    1381             CASE('slphynon') 
    1382 #if defined key_hadocc 
    1383                CALL ctl_stop( ' Trying to run slphynon observation operator', & 
    1384                   &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
    1385 #elif defined key_medusa 
    1386                ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
    1387                zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 
    1388 #elif defined key_fabm 
    1389                ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 
    1390                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
    1391                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
    1392 #else 
    1393                CALL ctl_stop( ' Trying to run slphynon observation operator', & 
    1394                   &           ' but no biogeochemical model appears to have been defined' ) 
    1395 #endif 
    1396                llog10 = .TRUE. 
    1397  
    1398             CASE('sspm') 
    1399 #if defined key_spm 
    1400                zsurfvar(:,:) = 0.0 
    1401                DO jn = 1, jp_spm 
    1402                   zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
    1403                END DO 
    1404 #else 
    1405                CALL ctl_stop( ' Trying to run sspm observation operator', & 
    1406                   &           ' but no spm model appears to have been defined' ) 
    1407 #endif 
    1408  
    1409             CASE('sfco2') 
    1410 #if defined key_hadocc 
    1411                zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
    1412                IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
    1413                   & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
    1414                   zsurfvar(:,:) = obfillflt 
    1415                   zsurfmask(:,:) = 0 
    1416                   CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
    1417                      &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
    1418                ENDIF 
    1419 #elif defined key_medusa && defined key_roam 
    1420                zsurfvar(:,:) = f2_fco2w(:,:) 
    1421 #elif defined key_fabm 
    1422                ! First, get pCO2 from FABM 
    1423                pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
    1424                zsurfvar(:,:) = pco2_3d(:,:,1) 
    1425                ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
    1426                ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
    1427                ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
    1428                ! and 
    1429                ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
    1430                ! Marine Chemistry, 2: 203-215. 
    1431                ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
    1432                ! not explicitly included - atmospheric pressure is not necessarily available so this is 
    1433                ! the best assumption. 
    1434                ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
    1435                ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
    1436                ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
    1437                ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
    1438                zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
    1439                   &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
    1440                   &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
    1441                   &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
    1442                   &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
    1443                   &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
    1444 #else 
    1445                CALL ctl_stop( ' Trying to run sfco2 observation operator', & 
    1446                   &           ' but no biogeochemical model appears to have been defined' ) 
    1447 #endif 
    1448  
    1449             CASE('spco2') 
    1450 #if defined key_hadocc 
    1451                zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
    1452                IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
    1453                   & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
    1454                   zsurfvar(:,:) = obfillflt 
    1455                   zsurfmask(:,:) = 0 
    1456                   CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
    1457                      &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
    1458                ENDIF 
    1459 #elif defined key_medusa && defined key_roam 
    1460                zsurfvar(:,:) = f2_pco2w(:,:) 
    1461 #elif defined key_fabm 
    1462                pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
    1463                zsurfvar(:,:) = pco2_3d(:,:,1) 
    1464 #else 
    1465                CALL ctl_stop( ' Trying to run spco2 observation operator', & 
    1466                   &           ' but no biogeochemical model appears to have been defined' ) 
    1467 #endif 
    1468  
    1469             CASE DEFAULT 
    1470  
    1471                CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 
    1472  
    1473             END SELECT 
    1474              
    1475             IF ( llog10 ) THEN 
    1476                ! Take the log10 where we can, otherwise exclude 
    1477                tiny = 1.0e-20 
    1478                WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
    1479                   zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
    1480                ELSEWHERE 
    1481                   zsurfvar(:,:)  = obfillflt 
    1482                   zsurfmask(:,:) = 0 
    1483                END WHERE 
    1484             ENDIF 
    1485  
    1486             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1487                &               nit000, idaystp, zsurfvar, zsurfmask,    & 
    1488                &               n2dintsurf(jtype), llnightav(jtype),     & 
    1489                &               ravglamscl(jtype), ravgphiscl(jtype),     & 
    1490                &               lfpindegs(jtype) ) 
    1491  
    1492          END DO 
    1493  
    1494          CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
    1495          CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    1496 #if defined key_fabm 
    1497          CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 
    1498 #endif 
    1499  
    1500       ENDIF 
     1211      IF ( nproftypes > 0 ) THEN  
     1212  
     1213         DO jtype = 1, nproftypes  
     1214  
     1215            ! Allocate local work arrays  
     1216            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )  
     1217            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )  
     1218            CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam     )  
     1219            CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi     )  
     1220            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim )      
     1221                                
     1222            ! Defaults which might change  
     1223            DO jvar = 1, profdataqc(jtype)%nvar  
     1224               zprofmask(:,:,:,jvar) = tmask(:,:,:)  
     1225               zglam(:,:,jvar)       = glamt(:,:)  
     1226               zgphi(:,:,jvar)       = gphit(:,:)  
     1227               zprofclim(:,:,:,jvar) = 0._wp  
     1228            END DO  
     1229  
     1230            SELECT CASE ( TRIM(cobstypesprof(jtype)) )  
     1231  
     1232            CASE('prof')  
     1233               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)  
     1234               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)  
     1235               IF ( ln_output_clim ) THEN            
     1236                  zprofclim(:,:,:,1) = tclim(:,:,:)  
     1237                  zprofclim(:,:,:,2) = sclim(:,:,:)  
     1238               ENDIF  
     1239                 
     1240            CASE('vel')  
     1241               zprofvar(:,:,:,1) = un(:,:,:)  
     1242               zprofvar(:,:,:,2) = vn(:,:,:)  
     1243               zprofmask(:,:,:,1) = umask(:,:,:)  
     1244               zprofmask(:,:,:,2) = vmask(:,:,:)  
     1245               zglam(:,:,1) = glamu(:,:)  
     1246               zglam(:,:,2) = glamv(:,:)  
     1247               zgphi(:,:,1) = gphiu(:,:)  
     1248               zgphi(:,:,2) = gphiv(:,:)  
     1249  
     1250            CASE('plchltot')  
     1251#if defined key_hadocc  
     1252               ! Chlorophyll from HadOCC  
     1253               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)  
     1254#elif defined key_medusa  
     1255               ! Add non-diatom and diatom chlorophyll from MEDUSA  
     1256               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)  
     1257#elif defined key_fabm  
     1258               ! Add all chlorophyll groups from ERSEM  
     1259               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &  
     1260                  & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)  
     1261#else  
     1262               CALL ctl_stop( ' Trying to run plchltot observation operator', &  
     1263                  &           ' but no biogeochemical model appears to have been defined' )  
     1264#endif  
     1265               ! Take the log10 where we can, otherwise exclude  
     1266               tiny = 1.0e-20  
     1267               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )  
     1268                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))  
     1269               ELSEWHERE  
     1270                  zprofvar(:,:,:,:)  = obfillflt  
     1271                  zprofmask(:,:,:,:) = 0  
     1272               END WHERE  
     1273               ! Mask out model below any excluded values,  
     1274               ! to avoid interpolation issues  
     1275               DO jvar = 1, profdataqc(jtype)%nvar  
     1276                 DO jj = 1, jpj  
     1277                    DO ji = 1, jpi  
     1278                       depth_loop: DO jk = 1, jpk  
     1279                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN  
     1280                             zprofmask(ji,jj,jk:jpk,jvar) = 0  
     1281                             EXIT depth_loop  
     1282                          ENDIF  
     1283                       END DO depth_loop  
     1284                    END DO  
     1285                 END DO  
     1286              END DO  
     1287  
     1288            CASE('pchltot')  
     1289#if defined key_hadocc  
     1290               ! Chlorophyll from HadOCC  
     1291               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)  
     1292#elif defined key_medusa  
     1293               ! Add non-diatom and diatom chlorophyll from MEDUSA  
     1294               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)  
     1295#elif defined key_fabm  
     1296               ! Add all chlorophyll groups from ERSEM  
     1297               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &  
     1298                  & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)  
     1299#else  
     1300               CALL ctl_stop( ' Trying to run pchltot observation operator', &  
     1301                  &           ' but no biogeochemical model appears to have been defined' )  
     1302#endif  
     1303  
     1304            CASE('pno3')  
     1305#if defined key_hadocc  
     1306               ! Dissolved inorganic nitrogen from HadOCC  
     1307               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)  
     1308#elif defined key_medusa  
     1309               ! Dissolved inorganic nitrogen from MEDUSA  
     1310               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)  
     1311#elif defined key_fabm  
     1312               ! Nitrate from ERSEM  
     1313               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)  
     1314#else  
     1315               CALL ctl_stop( ' Trying to run pno3 observation operator', &  
     1316                  &           ' but no biogeochemical model appears to have been defined' )  
     1317#endif  
     1318  
     1319            CASE('psi4')  
     1320#if defined key_hadocc  
     1321               CALL ctl_stop( ' Trying to run psi4 observation operator', &  
     1322                  &           ' but HadOCC does not simulate silicate' )  
     1323#elif defined key_medusa  
     1324               ! Silicate from MEDUSA  
     1325               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)  
     1326#elif defined key_fabm  
     1327               ! Silicate from ERSEM  
     1328               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)  
     1329#else  
     1330               CALL ctl_stop( ' Trying to run psi4 observation operator', &  
     1331                  &           ' but no biogeochemical model appears to have been defined' )  
     1332#endif  
     1333  
     1334            CASE('ppo4')  
     1335#if defined key_hadocc  
     1336               CALL ctl_stop( ' Trying to run ppo4 observation operator', &  
     1337                  &           ' but HadOCC does not simulate phosphate' )  
     1338#elif defined key_medusa  
     1339               CALL ctl_stop( ' Trying to run ppo4 observation operator', &  
     1340                  &           ' but MEDUSA does not simulate phosphate' )  
     1341#elif defined key_fabm  
     1342               ! Phosphate from ERSEM  
     1343               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)  
     1344#else  
     1345               CALL ctl_stop( ' Trying to run ppo4 observation operator', &  
     1346                  &           ' but no biogeochemical model appears to have been defined' )  
     1347#endif  
     1348  
     1349            CASE('pdic')  
     1350#if defined key_hadocc  
     1351               ! Dissolved inorganic carbon from HadOCC  
     1352               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)  
     1353#elif defined key_medusa  
     1354               ! Dissolved inorganic carbon from MEDUSA  
     1355               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)  
     1356#elif defined key_fabm  
     1357               ! Dissolved inorganic carbon from ERSEM  
     1358               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)  
     1359#else  
     1360               CALL ctl_stop( ' Trying to run pdic observation operator', &  
     1361                  &           ' but no biogeochemical model appears to have been defined' )  
     1362#endif  
     1363  
     1364            CASE('palk')  
     1365#if defined key_hadocc  
     1366               ! Alkalinity from HadOCC  
     1367               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)  
     1368#elif defined key_medusa  
     1369               ! Alkalinity from MEDUSA  
     1370               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)  
     1371#elif defined key_fabm  
     1372               ! Alkalinity from ERSEM  
     1373               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ta)  
     1374#else  
     1375               CALL ctl_stop( ' Trying to run palk observation operator', &  
     1376                  &           ' but no biogeochemical model appears to have been defined' )  
     1377#endif  
     1378  
     1379            CASE('pph')  
     1380#if defined key_hadocc  
     1381               CALL ctl_stop( ' Trying to run pph observation operator', &  
     1382                  &           ' but HadOCC has no pH diagnostic defined' )  
     1383#elif defined key_medusa && defined key_roam  
     1384               ! pH from MEDUSA  
     1385               zprofvar(:,:,:,1) = f3_pH(:,:,:)  
     1386#elif defined key_fabm  
     1387               ! pH from ERSEM  
     1388               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ph)  
     1389#else  
     1390               CALL ctl_stop( ' Trying to run pph observation operator', &  
     1391                  &           ' but no biogeochemical model appears to have been defined' )  
     1392#endif  
     1393  
     1394            CASE('po2')  
     1395#if defined key_hadocc  
     1396               CALL ctl_stop( ' Trying to run po2 observation operator', &  
     1397                  &           ' but HadOCC does not simulate oxygen' )  
     1398#elif defined key_medusa  
     1399               ! Oxygen from MEDUSA  
     1400               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)  
     1401#elif defined key_fabm  
     1402               ! Oxygen from ERSEM  
     1403               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)  
     1404#else  
     1405               CALL ctl_stop( ' Trying to run po2 observation operator', &  
     1406                  &           ' but no biogeochemical model appears to have been defined' )  
     1407#endif  
     1408  
     1409            CASE DEFAULT  
     1410               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )  
     1411  
     1412            END SELECT  
     1413  
     1414            DO jvar = 1, profdataqc(jtype)%nvar  
     1415               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &  
     1416                  &               nit000, idaystp, jvar,                   &  
     1417                  &               zprofvar(:,:,:,jvar),                    &  
     1418                  &               zprofclim(:,:,:,jvar),                   &  
     1419                  &               fsdept(:,:,:), fsdepw(:,:,:),            &   
     1420                  &               zprofmask(:,:,:,jvar),                   &  
     1421                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &  
     1422                  &               nn_1dint, nn_2dint_default,              &  
     1423                  &               kdailyavtypes = nn_profdavtypes )  
     1424            END DO  
     1425  
     1426            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )  
     1427            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )  
     1428            CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam     )  
     1429            CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi     )  
     1430            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim  )              
     1431  
     1432         END DO  
     1433  
     1434      ENDIF  
     1435  
     1436      IF ( nsurftypes > 0 ) THEN  
     1437           
     1438#if defined key_fabm  
     1439         CALL wrk_alloc( jpi, jpj, jpk, fabm_3d )  
     1440#endif  
     1441  
     1442         DO jtype = 1, nsurftypes  
     1443  
     1444            !Allocate local work arrays  
     1445            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  )  
     1446            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )           
     1447            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask )  
     1448            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     )  
     1449            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     )  
     1450  
     1451            !Defaults which might be changed  
     1452            DO jvar = 1, surfdataqc(jtype)%nvar              
     1453               zsurfmask(:,:,jvar) = tmask(:,:,1)  
     1454               zsurfclim(:,:,jvar) = 0._wp  
     1455               zglam(:,:,jvar) = glamt(:,:)  
     1456               zgphi(:,:,jvar) = gphit(:,:)     
     1457            END DO               
     1458            llog10 = .FALSE.  
     1459  
     1460            SELECT CASE ( TRIM(cobstypessurf(jtype)) )  
     1461            CASE('sst')  
     1462               zsurfvar(:,:,1) = tsn(:,:,1,jp_tem)  
     1463               IF ( ln_output_clim ) zsurfclim(:,:,1) = tclim(:,:,1)  
     1464            CASE('sla')  
     1465               zsurfvar(:,:,1) = sshn(:,:)  
     1466            CASE('sss')  
     1467               zsurfvar(:,:,1) = tsn(:,:,1,jp_sal)  
     1468               IF ( ln_output_clim ) zsurfclim(:,:,1) = sclim(:,:,1)                
     1469            CASE('ssv')  
     1470               zsurfvar(:,:,1) = un(:,:,1)  
     1471               zsurfvar(:,:,2) = vn(:,:,1)  
     1472               zsurfmask(:,:,1) = umask(:,:,1)  
     1473               zsurfmask(:,:,2) = vmask(:,:,1)     
     1474               zglam(:,:,1) = glamu(:,:)  
     1475               zglam(:,:,2) = glamv(:,:)                 
     1476               zgphi(:,:,1) = gphiu(:,:)  
     1477               zgphi(:,:,2) = gphiv(:,:)                    
     1478            CASE('sic')  
     1479               IF ( kstp == 0 ) THEN  
     1480                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN  
     1481                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &  
     1482                        &           'time-step but some obs are valid then.' )  
     1483                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &  
     1484                        &           ' sea-ice concentration obs will be missed'  
     1485                  ENDIF  
     1486               ELSE  
     1487#if defined key_cice  
     1488                  zsurfvar(:,:,1) = fr_i(:,:)  
     1489#elif defined key_lim2 || defined key_lim3  
     1490                  zsurfvar(:,:,1) = 1._wp - frld(:,:)  
     1491#else  
     1492               CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', &  
     1493                  &           ' but no sea-ice model appears to have been defined' )  
     1494#endif  
     1495               ENDIF  
     1496  
     1497            CASE('sit')  
     1498               IF ( kstp == 0 ) THEN  
     1499                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN  
     1500                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &  
     1501                        &           'time-step but some obs are valid then.' )  
     1502                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &  
     1503                        &           ' sea-ice thickness obs will be missed and QC flag set to 4'  
     1504                  ENDIF                                   
     1505               ELSE         
     1506#if defined key_cice  
     1507                  zsurfvar(:,:,1) = thick_i(:,:)  
     1508#elif defined key_lim2 || defined key_lim3  
     1509                  CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' )  
     1510#else  
     1511                  CALL ctl_stop( ' Trying to run sea-ice thickness observation operator', &  
     1512                     &           ' but no sea-ice model appears to have been defined' )  
     1513#endif  
     1514               ENDIF  
     1515  
     1516            CASE('fbd')  
     1517               IF ( kstp == 0 ) THEN  
     1518                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN  
     1519                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &  
     1520                        &           'time-step but some obs are valid then.' )  
     1521                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &  
     1522                        &           ' sea-ice freeboard obs will be missed and QC flag set to 4'  
     1523                  ENDIF                                   
     1524               ELSE         
     1525#if defined key_cice  
     1526                  zsurfvar(:,:) = thick_i(:,:)  
     1527#elif defined key_lim2 || defined key_lim3  
     1528                  CALL ctl_stop( ' No sea-ice freeboard observation operator defined for LIM model' )  
     1529#else  
     1530                  CALL ctl_stop( ' Trying to run sea-ice freeboard observation operator', &  
     1531                     &           ' but no sea-ice model appears to have been defined' )  
     1532#endif  
     1533               ENDIF  
     1534  
     1535            CASE('slchltot')  
     1536#if defined key_hadocc  
     1537               ! Surface chlorophyll from HadOCC  
     1538               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1)  
     1539#elif defined key_medusa  
     1540               ! Add non-diatom and diatom surface chlorophyll from MEDUSA  
     1541               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)  
     1542#elif defined key_fabm  
     1543               ! Add all surface chlorophyll groups from ERSEM  
     1544               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &  
     1545                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)  
     1546#else  
     1547               CALL ctl_stop( ' Trying to run slchltot observation operator', &  
     1548                  &           ' but no biogeochemical model appears to have been defined' )  
     1549#endif  
     1550               llog10 = .TRUE.  
     1551  
     1552            CASE('slchldia')  
     1553#if defined key_hadocc  
     1554               CALL ctl_stop( ' Trying to run slchldia observation operator', &  
     1555                  &           ' but HadOCC does not explicitly simulate diatoms' )  
     1556#elif defined key_medusa  
     1557               ! Diatom surface chlorophyll from MEDUSA  
     1558               zsurfvar(:,:,1) = trn(:,:,1,jpchd)  
     1559#elif defined key_fabm  
     1560               ! Diatom surface chlorophyll from ERSEM  
     1561               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)  
     1562#else  
     1563               CALL ctl_stop( ' Trying to run slchldia observation operator', &  
     1564                  &           ' but no biogeochemical model appears to have been defined' )  
     1565#endif  
     1566               llog10 = .TRUE.  
     1567  
     1568            CASE('slchlnon')  
     1569#if defined key_hadocc  
     1570               CALL ctl_stop( ' Trying to run slchlnon observation operator', &  
     1571                  &           ' but HadOCC does not explicitly simulate non-diatoms' )  
     1572#elif defined key_medusa  
     1573               ! Non-diatom surface chlorophyll from MEDUSA  
     1574               zsurfvar(:,:,1) = trn(:,:,1,jpchn)  
     1575#elif defined key_fabm  
     1576               ! Add all non-diatom surface chlorophyll groups from ERSEM  
     1577               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &  
     1578                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)  
     1579#else  
     1580               CALL ctl_stop( ' Trying to run slchlnon observation operator', &  
     1581                  &           ' but no biogeochemical model appears to have been defined' )  
     1582#endif  
     1583               llog10 = .TRUE.  
     1584  
     1585            CASE('slchldin')  
     1586#if defined key_hadocc  
     1587               CALL ctl_stop( ' Trying to run slchldin observation operator', &  
     1588                  &           ' but HadOCC does not explicitly simulate dinoflagellates' )  
     1589#elif defined key_medusa  
     1590               CALL ctl_stop( ' Trying to run slchldin observation operator', &  
     1591                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' )  
     1592#elif defined key_fabm  
     1593               ! Dinoflagellate surface chlorophyll from ERSEM  
     1594               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)  
     1595#else  
     1596               CALL ctl_stop( ' Trying to run slchldin observation operator', &  
     1597                  &           ' but no biogeochemical model appears to have been defined' )  
     1598#endif  
     1599               llog10 = .TRUE.  
     1600  
     1601            CASE('slchlmic')  
     1602#if defined key_hadocc  
     1603               CALL ctl_stop( ' Trying to run slchlmic observation operator', &  
     1604                  &           ' but HadOCC does not explicitly simulate microphytoplankton' )  
     1605#elif defined key_medusa  
     1606               CALL ctl_stop( ' Trying to run slchlmic observation operator', &  
     1607                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' )  
     1608#elif defined key_fabm  
     1609               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM  
     1610               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)  
     1611#else  
     1612               CALL ctl_stop( ' Trying to run slchlmic observation operator', &  
     1613                  &           ' but no biogeochemical model appears to have been defined' )  
     1614#endif  
     1615               llog10 = .TRUE.  
     1616  
     1617            CASE('slchlnan')  
     1618#if defined key_hadocc  
     1619               CALL ctl_stop( ' Trying to run slchlnan observation operator', &  
     1620                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' )  
     1621#elif defined key_medusa  
     1622               CALL ctl_stop( ' Trying to run slchlnan observation operator', &  
     1623                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' )  
     1624#elif defined key_fabm  
     1625               ! Nanophytoplankton surface chlorophyll from ERSEM  
     1626               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)  
     1627#else  
     1628               CALL ctl_stop( ' Trying to run slchlnan observation operator', &  
     1629                  &           ' but no biogeochemical model appears to have been defined' )  
     1630#endif  
     1631               llog10 = .TRUE.  
     1632  
     1633            CASE('slchlpic')  
     1634#if defined key_hadocc  
     1635               CALL ctl_stop( ' Trying to run slchlpic observation operator', &  
     1636                  &           ' but HadOCC does not explicitly simulate picophytoplankton' )  
     1637#elif defined key_medusa  
     1638               CALL ctl_stop( ' Trying to run slchlpic observation operator', &  
     1639                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' )  
     1640#elif defined key_fabm  
     1641               ! Picophytoplankton surface chlorophyll from ERSEM  
     1642               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)  
     1643#else  
     1644               CALL ctl_stop( ' Trying to run slchlpic observation operator', &  
     1645                  &           ' but no biogeochemical model appears to have been defined' )  
     1646#endif  
     1647               llog10 = .TRUE.  
     1648  
     1649            CASE('schltot')  
     1650#if defined key_hadocc  
     1651               ! Surface chlorophyll from HadOCC  
     1652               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1)  
     1653#elif defined key_medusa  
     1654               ! Add non-diatom and diatom surface chlorophyll from MEDUSA  
     1655               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)  
     1656#elif defined key_fabm  
     1657               ! Add all surface chlorophyll groups from ERSEM  
     1658               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &  
     1659                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)  
     1660#else  
     1661               CALL ctl_stop( ' Trying to run schltot observation operator', &  
     1662                  &           ' but no biogeochemical model appears to have been defined' )  
     1663#endif  
     1664  
     1665            CASE('slphytot')  
     1666#if defined key_hadocc  
     1667               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio  
     1668               zsurfvar(:,:,1) = trn(:,:,1,jp_had_phy) * c2n_p  
     1669#elif defined key_medusa  
     1670               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA  
     1671               ! multiplied by C:N ratio for each  
     1672               zsurfvar(:,:,1) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)  
     1673#elif defined key_fabm  
     1674               ! Add all surface phytoplankton carbon groups from ERSEM  
     1675               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &  
     1676                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)  
     1677#else  
     1678               CALL ctl_stop( ' Trying to run slphytot observation operator', &  
     1679                  &           ' but no biogeochemical model appears to have been defined' )  
     1680#endif  
     1681               llog10 = .TRUE.  
     1682  
     1683            CASE('slphydia')  
     1684#if defined key_hadocc  
     1685               CALL ctl_stop( ' Trying to run slphydia observation operator', &  
     1686                  &           ' but HadOCC does not explicitly simulate diatoms' )  
     1687#elif defined key_medusa  
     1688               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio  
     1689               zsurfvar(:,:,1) = trn(:,:,1,jpphd) * xthetapd  
     1690#elif defined key_fabm  
     1691               ! Diatom surface phytoplankton carbon from ERSEM  
     1692               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)  
     1693#else  
     1694               CALL ctl_stop( ' Trying to run slphydia observation operator', &  
     1695                  &           ' but no biogeochemical model appears to have been defined' )  
     1696#endif  
     1697               llog10 = .TRUE.  
     1698  
     1699            CASE('slphynon')  
     1700#if defined key_hadocc  
     1701               CALL ctl_stop( ' Trying to run slphynon observation operator', &  
     1702                  &           ' but HadOCC does not explicitly simulate non-diatoms' )  
     1703#elif defined key_medusa  
     1704               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio  
     1705               zsurfvar(:,:,1) = trn(:,:,1,jpphn) * xthetapn  
     1706#elif defined key_fabm  
     1707               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM  
     1708               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &  
     1709                  & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)  
     1710#else  
     1711               CALL ctl_stop( ' Trying to run slphynon observation operator', &  
     1712                  &           ' but no biogeochemical model appears to have been defined' )  
     1713#endif  
     1714               llog10 = .TRUE.  
     1715  
     1716            CASE('sspm')  
     1717#if defined key_spm  
     1718               zsurfvar(:,:,1) = 0.0  
     1719               DO jn = 1, jp_spm  
     1720                  zsurfvar(:,:,1) = zsurfvar(:,:,1) + trn(:,:,1,jn)   ! sum SPM sizes  
     1721               END DO  
     1722#else  
     1723               CALL ctl_stop( ' Trying to run sspm observation operator', &  
     1724                  &           ' but no spm model appears to have been defined' )  
     1725#endif  
     1726  
     1727            CASE('skd490')  
     1728#if defined key_hadocc  
     1729               CALL ctl_stop( ' Trying to run skd490 observation operator', &  
     1730                  &           ' but HadOCC does not explicitly simulate Kd490' )  
     1731#elif defined key_medusa  
     1732               CALL ctl_stop( ' Trying to run skd490 observation operator', &  
     1733                  &           ' but MEDUSA does not explicitly simulate Kd490' )  
     1734#elif defined key_fabm  
     1735               ! light_Kd_band3 diagnostic variable if using spectral optical model  
     1736               ! light_xEPS diagnostic variable if using standard ERSEM light model  
     1737               IF ( jp_fabm_kd490 /= -1 ) THEN  
     1738                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_kd490)  
     1739               ELSEIF ( jp_fabm_xeps /= -1 ) THEN  
     1740                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_xeps)  
     1741               ELSE  
     1742                  CALL ctl_stop( ' Trying to run skd490 observation operator', &  
     1743                     &           ' but cannot access Kd490 from ERSEM' )  
     1744               ENDIF  
     1745               zsurfvar(:,:,1) = fabm_3d(:,:,1)  
     1746#else  
     1747               CALL ctl_stop( ' Trying to run skd490 observation operator', &  
     1748                  &           ' but no biogeochemical model appears to have been defined' )  
     1749#endif  
     1750  
     1751            CASE('sfco2')  
     1752#if defined key_hadocc  
     1753               zsurfvar(:,:,1) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC  
     1754               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &  
     1755                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN  
     1756                  zsurfvar(:,:) = obfillflt  
     1757                  zsurfmask(:,:,1) = 0  
     1758                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &  
     1759                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )  
     1760               ENDIF  
     1761#elif defined key_medusa && defined key_roam  
     1762               zsurfvar(:,:,1) = f2_fco2w(:,:)  
     1763#elif defined key_fabm  
     1764               ! First, get pCO2 from FABM  
     1765               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc)  
     1766               zsurfvar(:,:,1) = fabm_3d(:,:,1)  
     1767               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:  
     1768               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems  
     1769               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.  
     1770               ! and  
     1771               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,  
     1772               ! Marine Chemistry, 2: 203-215.  
     1773               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so  
     1774               ! not explicitly included - atmospheric pressure is not necessarily available so this is  
     1775               ! the best assumption.  
     1776               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice  
     1777               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)  
     1778               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal  
     1779               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.  
     1780               zsurfvar(:,:,1) = zsurfvar(:,:,1) * EXP((-1636.75 + &  
     1781                  &              12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &  
     1782                  &              0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &  
     1783                  &              0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &  
     1784                  &              2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &  
     1785                  &              (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))  
     1786#else  
     1787               CALL ctl_stop( ' Trying to run sfco2 observation operator', &  
     1788                  &           ' but no biogeochemical model appears to have been defined' )  
     1789#endif  
     1790  
     1791            CASE('spco2')  
     1792#if defined key_hadocc  
     1793               zsurfvar(:,:,1) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC  
     1794               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &  
     1795                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN  
     1796                  zsurfvar(:,:,1) = obfillflt  
     1797                  zsurfmask(:,:,1) = 0  
     1798                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &  
     1799                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )  
     1800               ENDIF  
     1801#elif defined key_medusa && defined key_roam  
     1802               zsurfvar(:,:,1) = f2_pco2w(:,:)  
     1803#elif defined key_fabm  
     1804               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc)  
     1805               zsurfvar(:,:,1) = fabm_3d(:,:,1)  
     1806#else  
     1807               CALL ctl_stop( ' Trying to run spco2 observation operator', &  
     1808                  &           ' but no biogeochemical model appears to have been defined' )  
     1809#endif  
     1810  
     1811            CASE DEFAULT  
     1812  
     1813               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )  
     1814  
     1815            END SELECT  
     1816              
     1817            IF ( llog10 ) THEN  
     1818               ! Take the log10 where we can, otherwise exclude  
     1819               tiny = 1.0e-20  
     1820               WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt )  
     1821                  zsurfvar(:,:,1)  = LOG10(zsurfvar(:,:,1))  
     1822               ELSEWHERE  
     1823                  zsurfvar(:,:,1)  = obfillflt  
     1824                  zsurfmask(:,:,1) = 0  
     1825               END WHERE  
     1826            ENDIF  
     1827  
     1828            DO jvar = 1, surfdataqc(jtype)%nvar  
     1829  
     1830               IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND. &  
     1831                     &  ln_time_mean_sla_bkg ) THEN  
     1832                  !Number of time-steps in meaning period  
     1833                  imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt )  
     1834                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &  
     1835                     &               nit000, idaystp, jvar,                   &  
     1836                     &               zsurfvar(:,:,jvar),                      &  
     1837                     &               zsurfclim(:,:,jvar),                     &  
     1838                     &               zsurfmask(:,:,jvar),                     &  
     1839                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                       
     1840                     &               n2dintsurf(jtype), llnightav(jtype),     &  
     1841                     &               ravglamscl(jtype), ravgphiscl(jtype),    &  
     1842                     &               lfpindegs(jtype), kmeanstp = imeanstp )  
     1843  
     1844               ELSE  
     1845                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &  
     1846                     &               nit000, idaystp, jvar,                   &  
     1847                     &               zsurfvar(:,:,jvar),                      &  
     1848                     &               zsurfclim(:,:,jvar),                     &  
     1849                     &               zsurfmask(:,:,jvar),                     &  
     1850                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                                            
     1851                     &               n2dintsurf(jtype), llnightav(jtype),     &  
     1852                     &               ravglamscl(jtype), ravgphiscl(jtype),    &  
     1853                     &               lfpindegs(jtype) )  
     1854               ENDIF  
     1855  
     1856            END DO  
     1857  
     1858            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  )  
     1859            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )                    
     1860            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask )  
     1861            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     )  
     1862            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     )           
     1863  
     1864         END DO  
     1865#if defined key_fabm  
     1866         CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d )  
     1867#endif  
     1868  
     1869      ENDIF  
    15011870 
    15021871   END SUBROUTINE dia_obs 
    … …  
    15491918                  & ) 
    15501919 
    1551                CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
     1920               CALL obs_rotvel_pro( profdataqc(jtype), nn_2dint_default, zu, zv ) 
    15521921 
    15531922               DO jo = 1, profdataqc(jtype)%nprof 
    … …  
    15821951 
    15831952         DO jtype = 1, nsurftypes 
     1953 
     1954            IF ( TRIM(cobstypessurf(jtype)) == 'vel' ) THEN  
     1955               ! For velocity data, rotate the model velocities to N/S, E/W  
     1956               ! using the compressed data structure.  
     1957               ALLOCATE( &  
     1958                  & zu(surfdataqc(jtype)%nsurf), &  
     1959                  & zv(surfdataqc(jtype)%nsurf)  &  
     1960                  & )  
     1961 
     1962               CALL obs_rotvel_surf( surfdataqc(jtype), nn_2dint_default, zu, zv )  
     1963 
     1964               DO jo = 1, surfdataqc(jtype)%nsurf  
     1965                  surfdataqc(jtype)%rmod(jo,1) = zu(jo)  
     1966                  surfdataqc(jtype)%rmod(jo,2) = zv(jo)  
     1967               END DO  
     1968 
     1969               DEALLOCATE( zu )  
     1970               DEALLOCATE( zv )  
     1971            END IF  
    15841972 
    15851973            CALL obs_surf_decompress( surfdataqc(jtype), & 

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

    r11203 r15764  
    1313      !! ** Method  : Straightforward search 
    1414      !! 
    15       !! ** Action  :  
     15      !! ** Action  : Will return level associated with T-point below the obs  
     16      !!              depth, except when observation is in the top box will   
     17      !!              return level 2. Also, if obs depth greater than depth   
     18      !!              of last wet T-point (kpk-1) will return level kpk. 
    1619      !! 
    1720      !! History : 
    … …  
    4346      DO ji = 1, kobs  
    4447         kobsk(ji) = 1 
    45          depk: DO jk = 2, kgrd 
    46             IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk 
     48         depk: DO jk = 2, kgrd-1  
     49            IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk 
    4750         END DO depk 
    4851         kobsk(ji) = jk 

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

    r11203 r15764  
    3737   USE obs_grid,      ONLY : &  
    3838      & obs_level_search      
     39#if defined key_cice  
     40   USE ice_constants, ONLY : &    ! For conversion from sea ice freeboard to thickness  
     41      & rhos, rhoi, rhow  
     42#endif 
    3943 
    4044   IMPLICIT NONE 
    … …  
    6064 
    6165 
    62    SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
    63       &                     kit000, kdaystp, kvar,       & 
    64       &                     pvar, pgdept, pgdepw,        & 
    65       &                     pmask,                       &   
    66       &                     plam, pphi,                  & 
     66   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &  
     67      &                     kit000, kdaystp, kvar,       &  
     68      &                     pvar, pclim,                 &  
     69      &                     pgdept, pgdepw, pmask,       &    
     70      &                     plam, pphi,                  &  
    6771      &                     k1dint, k2dint, kdailyavtypes ) 
    6872 
    … …  
    137141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    138142         & pvar,   &                 ! Model field for variable 
     143         & pclim,  &                 ! Climatology field for variable 
    139144         & pmask                     ! Land-sea mask for variable 
    140145      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    … …  
    172177      REAL(KIND=wp), DIMENSION(kpk) :: & 
    173178         & zobsk,    & 
    174          & zobs2k 
     179         & zobs2k,   &  
     180         & zclm2k 
    175181      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    176182         & zweig1, & 
    … …  
    178184      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    179185         & zmask,  & 
     186         & zclim,  & 
    180187         & zint,   & 
    181188         & zinm,   & 
    … …  
    187194      REAL(KIND=wp), DIMENSION(1) :: zmsk 
    188195      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
     196      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim 
    189197 
    190198      LOGICAL :: ld_dailyav 
    … …  
    262270         & ) 
    263271 
     272      IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) ) 
     273 
    264274      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    265275         iobs = jobs - prodatqc%nprofup 
    … …  
    285295      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )  
    286296      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )  
     297 
     298      IF ( prodatqc%lclim ) THEN  
     299         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )              
     300      ENDIF 
    287301 
    288302      ! At the end of the day also get interpolated means 
    … …  
    349363                  inum_obs = iend - ista + 1  
    350364                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     365                  IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 
    351366 
    352367                  DO iin=1,2  
    … …  
    358373                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    359374                              &     zmask(iin,ijn,:,iobs))  
     375 
     376                           IF ( prodatqc%lclim ) THEN  
     377                              CALL obs_int_z1d_spl( kpk, &   
     378                                 &     zclim(iin,ijn,:,iobs), &   
     379                                 &     zclm2k, zgdept(iin,ijn,:,iobs), &   
     380                                 &     zmask(iin,ijn,:,iobs))   
     381                           ENDIF 
    360382                        ENDIF  
    361383        
    … …  
    372394                           &    zmask(iin,ijn,:,iobs))  
    373395        
     396                        IF ( prodatqc%lclim ) THEN  
     397                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &   
     398                              &    prodatqc%var(kvar)%vdep(ista:iend), &   
     399                              &    zclim(iin,ijn,:,iobs), &   
     400                              &    zclm2k, interp_corner_clim(iin,ijn,:), &   
     401                              &    zgdept(iin,ijn,:,iobs), &   
     402                              &    zmask(iin,ijn,:,iobs))   
     403                        ENDIF 
    374404                     ENDDO  
    375405                  ENDDO  
    … …  
    386416               inum_obs = iend - ista + 1  
    387417               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     418               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 
    388419               DO iin=1,2   
    389420                  DO ijn=1,2  
    … …  
    395426                           &    zmask(iin,ijn,:,iobs))  
    396427   
     428                        IF ( prodatqc%lclim ) THEN  
     429                           CALL obs_int_z1d_spl( kpk, &   
     430                              &    zclim(iin,ijn,:,iobs),&   
     431                              &    zclm2k, zgdept(iin,ijn,:,iobs), &   
     432                              &    zmask(iin,ijn,:,iobs))   
     433                        ENDIF 
    397434                     ENDIF  
    398435        
    … …  
    408445                         &          zgdept(iin,ijn,:,iobs),         &  
    409446                         &          zmask(iin,ijn,:,iobs) )       
    410           
     447 
     448                     IF ( prodatqc%lclim ) THEN  
     449                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &   
     450                            &          prodatqc%var(kvar)%vdep(ista:iend),     &   
     451                            &          zclim(iin,ijn,:,iobs),  &   
     452                            &          zclm2k,interp_corner_clim(iin,ijn,:), &   
     453                            &          zgdept(iin,ijn,:,iobs), &   
     454                            &          zmask(iin,ijn,:,iobs) )     
     455                     ENDIF 
    411456                  ENDDO  
    412457               ENDDO  
    … …  
    451496                  &              prodatqc%var(kvar)%vmod(iend:iend) )  
    452497 
     498               IF ( prodatqc%lclim ) THEN  
     499                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), &   
     500                     &              prodatqc%var(kvar)%vclm(iend:iend) )  
     501               ENDIF 
     502 
    453503                  ! Set QC flag for any observations found below the bottom 
    454504                  ! needed as the check here is more strict than that in obs_prep 
    … …  
    458508  
    459509            DEALLOCATE(interp_corner,iv_indic)  
     510            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim ) 
    460511           
    461512         ENDIF 
    … …  
    475526         & ) 
    476527 
     528      IF ( prodatqc%lclim ) DEALLOCATE( zclim ) 
     529 
    477530      ! At the end of the day also get interpolated means 
    478531      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    … …  
    486539   END SUBROUTINE obs_prof_opt 
    487540 
    488    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
    489       &                     kit000, kdaystp, psurf, psurfmask,   & 
    490       &                     k2dint, ldnightav, plamscl, pphiscl, & 
    491       &                     lindegrees ) 
     541   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &  
     542      &                     kit000, kdaystp, kvar, &  
     543      &                     psurf, pclim, psurfmask, &  
     544      &                     plam, pphi, &  
     545      &                     k2dint, ldnightav, plamscl, pphiscl, &  
     546      &                     lindegrees, kmeanstp ) 
    492547 
    493548      !!----------------------------------------------------------------------- 
    … …  
    522577      !!      ! 15-02 (M. Martin) Combined routine for surface types 
    523578      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
     579      !!      ! 20-08 (M. Martin) Added surface velocity options 
    524580      !!----------------------------------------------------------------------- 
    525581 
    … …  
    538594                                       !   (kit000-1 = restart time) 
    539595      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
     596      INTEGER, INTENT(IN) :: kvar      ! Number of variables in surfdataqc 
    540597      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     598      INTEGER, INTENT(IN), OPTIONAL :: &  
     599                             kmeanstp  ! Number of time steps for the time meaning  
     600                                       ! Averaging is triggered if present and greater than one 
    541601      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    542602         & psurf,  &                   ! Model surface field 
     603         & pclim,  &                   ! Climatological surface field 
    543604         & psurfmask                   ! Land-sea mask 
     605      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     606         & plam,   &                   ! Model longitudes for variable  
     607         & pphi                        ! Model latitudes for variable 
    544608      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
    545609      REAL(KIND=wp), INTENT(IN) :: & 
    … …  
    559623      INTEGER :: imodi, imodj 
    560624      INTEGER :: idayend 
     625      INTEGER :: imeanend 
    561626      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    562627         & igrdi,   & 
    … …  
    569634      REAL(wp) :: zlam 
    570635      REAL(wp) :: zphi 
    571       REAL(wp), DIMENSION(1) :: zext, zobsmask 
     636      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 
    572637      REAL(wp) :: zdaystp 
     638      REAL(wp) :: zmeanstp 
    573639      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    574640         & zweig,  & 
    … …  
    577643         & zsurfm, & 
    578644         & zsurftmp, & 
     645         & zclim,  & 
    579646         & zglam,  & 
    580647         & zgphi,  & 
    … …  
    586653         & zouttmp, & 
    587654         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
     655 
     656      LOGICAL :: l_timemean 
    588657 
    589658      !------------------------------------------------------------------------ 
    … …  
    594663      isurf = surfdataqc%nsstp(inrc) 
    595664 
     665      l_timemean = .FALSE.  
     666      IF ( PRESENT( kmeanstp ) ) THEN  
     667         IF ( kmeanstp > 1 ) l_timemean = .TRUE.  
     668      ENDIF 
     669 
    596670      ! Work out the maximum footprint size for the  
    597671      ! interpolation/averaging in model grid-points - has to be even. 
    … …  
    599673      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
    600674 
     675      IF ( l_timemean ) THEN  
     676         ! Initialize time mean for first timestep  
     677         imeanend = MOD( kt - kit000 + 1, kmeanstp )  
     678         IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend  
     679 
     680         ! Added kt == 0 test to catch restart case   
     681         IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN  
     682            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt  
     683            DO jj = 1, jpj  
     684               DO ji = 1, jpi  
     685                  surfdataqc%vdmean(ji,jj,kvar) = 0.0  
     686               END DO  
     687            END DO  
     688         ENDIF  
     689 
     690         ! On each time-step, increment the field for computing time mean  
     691         IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt  
     692         DO jj = 1, jpj  
     693            DO ji = 1, jpi  
     694               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &  
     695                  &                            + psurf(ji,jj)  
     696            END DO  
     697         END DO  
     698 
     699         ! Compute the time mean at the end of time period  
     700         IF ( imeanend == 0 ) THEN  
     701            zmeanstp = 1.0 / REAL( kmeanstp )  
     702            IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp  
     703            DO jj = 1, jpj  
     704               DO ji = 1, jpi  
     705                  surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &  
     706                     &                            * zmeanstp  
     707               END DO  
     708            END DO  
     709         ENDIF  
     710      ENDIF !l_timemean 
    601711 
    602712      IF ( ldnightav ) THEN 
    … …  
    665775 
    666776      ALLOCATE( & 
    667          & zweig(imaxifp,imaxjfp,1),      & 
    668          & igrdi(imaxifp,imaxjfp,isurf), & 
    669          & igrdj(imaxifp,imaxjfp,isurf), & 
    670          & zglam(imaxifp,imaxjfp,isurf), & 
    671          & zgphi(imaxifp,imaxjfp,isurf), & 
    672          & zmask(imaxifp,imaxjfp,isurf), & 
    673          & zsurf(imaxifp,imaxjfp,isurf), & 
    674          & zsurftmp(imaxifp,imaxjfp,isurf),  & 
    675          & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
    676          & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
    677          & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
    678          & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    679          & ) 
     777         & zweig(imaxifp,imaxjfp,1),       &  
     778         & igrdi(imaxifp,imaxjfp,isurf),   &  
     779         & igrdj(imaxifp,imaxjfp,isurf),   &  
     780         & zglam(imaxifp,imaxjfp,isurf),   &  
     781         & zgphi(imaxifp,imaxjfp,isurf),   &  
     782         & zmask(imaxifp,imaxjfp,isurf),   &  
     783         & zsurf(imaxifp,imaxjfp,isurf),   &  
     784         & zsurftmp(imaxifp,imaxjfp,isurf) &  
     785         & )  
     786 
     787      IF ( k2dint > 4 ) THEN  
     788         ALLOCATE( &  
     789            & zglamf(imaxifp+1,imaxjfp+1,isurf),  &  
     790            & zgphif(imaxifp+1,imaxjfp+1,isurf),  &  
     791            & igrdip1(imaxifp+1,imaxjfp+1,isurf), &  
     792            & igrdjp1(imaxifp+1,imaxjfp+1,isurf)  &  
     793            & )  
     794      ENDIF  
     795       
     796      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )  
    680797 
    681798      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    … …  
    695812               IF ( imodj > jpjglo ) imodj = jpjglo 
    696813 
    697                igrdip1(ji+1,jj+1,iobs) = imodi 
    698                igrdjp1(ji+1,jj+1,iobs) = imodj 
     814               IF ( k2dint > 4 ) THEN  
     815                  igrdip1(ji+1,jj+1,iobs) = imodi  
     816                  igrdjp1(ji+1,jj+1,iobs) = imodj  
     817               ENDIF 
    699818                
    700819               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
    … …  
    713832      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    714833         &                  igrdi, igrdj, psurfmask, zmask ) 
    715       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    716          &                  igrdi, igrdj, psurf, zsurf ) 
    717       CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
    718          &                  igrdip1, igrdjp1, glamf, zglamf ) 
    719       CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
    720          &                  igrdip1, igrdjp1, gphif, zgphif ) 
     834 
     835      ! At the end of the averaging period get interpolated means  
     836      IF ( l_timemean ) THEN  
     837         IF ( imeanend == 0 ) THEN  
     838            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )  
     839            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt  
     840            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &  
     841               &                  igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )  
     842         ENDIF  
     843      ELSE  
     844         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &  
     845            &                  igrdi, igrdj, psurf, zsurf )  
     846      ENDIF  
     847 
     848      IF ( k2dint > 4 ) THEN           
     849         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &  
     850            &                  igrdip1, igrdjp1, glamf, zglamf )  
     851         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &  
     852            &                  igrdip1, igrdjp1, gphif, zgphif )  
     853      ENDIF 
     854 
     855      IF ( surfdataqc%lclim ) THEN   
     856         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &  
     857            &                  igrdi, igrdj, pclim, zclim )  
     858      ENDIF 
    721859 
    722860      ! At the end of the day get interpolated means 
    … …  
    758896         zphi = surfdataqc%rphi(jobs) 
    759897 
    760          IF ( ldnightav .AND. idayend == 0 ) THEN 
    761             ! Night-time averaged data 
    762             zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    763          ELSE 
    764             zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     898         IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN  
     899            ! Night-time or N=kmeanstp timestep averaged data  
     900            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)  
     901         ELSE  
     902            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)  
     903         ENDIF  
     904 
     905         IF (   ( .NOT. l_timemean ) .OR. &   
     906             &  ( l_timemean .AND. imeanend == 0) ) THEN  
     907            IF ( k2dint <= 4 ) THEN  
     908 
     909               ! Get weights to interpolate the model value to the observation point  
     910               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &  
     911                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
     912                  &                   zmask(:,:,iobs), zweig, zobsmask )  
     913 
     914               ! Interpolate the model value to the observation point   
     915               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )  
     916 
     917               IF ( surfdataqc%lclim ) THEN    
     918                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )  
     919               ENDIF  
     920 
     921 
     922            ELSE  
     923 
     924               ! Get weights to average the model SLA to the observation footprint  
     925               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &  
     926                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
     927                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &  
     928                  &                   zmask(:,:,iobs), plamscl, pphiscl, &  
     929                  &                   lindegrees, zweig, zobsmask )  
     930 
     931               ! Average the model SST to the observation footprint  
     932               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &  
     933                  &              zweig, zsurftmp(:,:,iobs), zext )  
     934 
     935               IF ( surfdataqc%lclim ) THEN    
     936                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &  
     937                     &              zweig, zclim(:,:,iobs), zclm )  
     938               ENDIF  
     939 
     940            ENDIF  
     941 
     942            IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN  
     943               ! ... Remove the MDT from the SSH at the observation point to get the SLA  
     944               surfdataqc%rext(jobs,1) = zext(1)  
     945               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)  
     946            ELSE  
     947               surfdataqc%rmod(jobs,kvar) = zext(1)  
     948            ENDIF  
     949 
     950            IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)  
     951 
     952            IF ( zext(1) == obfillflt ) THEN  
     953               ! If the observation value is a fill value, set QC flag to bad  
     954               surfdataqc%nqc(jobs) = 4  
     955            ENDIF           
     956         ENDIF  
     957 
     958         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)  
     959  
     960#if defined key_cice  
     961         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN  
     962            ! Convert radar freeboard to true freeboard (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum compared to snow (3.0e8 vs 2.4e8 m/s))  
     963            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)   
     964            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)  
     965            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad  
     966            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN  
     967               surfdataqc%nqc(jobs) = 4  
     968            ENDIF             
     969            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016)  
     970            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ &  
     971                                      (rhow - rhoi)  
    765972         ENDIF 
    766973 
    767          IF ( k2dint <= 4 ) THEN 
    768  
    769             ! Get weights to interpolate the model value to the observation point 
    770             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    771                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    772                &                   zmask(:,:,iobs), zweig, zobsmask ) 
    773  
    774             ! Interpolate the model value to the observation point  
    775             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
    776  
    777          ELSE 
    778  
    779             ! Get weights to average the model SLA to the observation footprint 
    780             CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
    781                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    782                &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
    783                &                   zmask(:,:,iobs), plamscl, pphiscl, & 
    784                &                   lindegrees, zweig, zobsmask ) 
    785  
    786             ! Average the model SST to the observation footprint 
    787             CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
    788                &              zweig, zsurftmp(:,:,iobs),  zext ) 
    789  
    790          ENDIF 
    791  
    792          IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
    793             ! ... Remove the MDT from the SSH at the observation point to get the SLA 
    794             surfdataqc%rext(jobs,1) = zext(1) 
    795             surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
    796          ELSE 
    797             surfdataqc%rmod(jobs,1) = zext(1) 
    798          ENDIF 
    799           
    800          IF ( zext(1) == obfillflt ) THEN 
    801             ! If the observation value is a fill value, set QC flag to bad 
    802             surfdataqc%nqc(jobs) = 4 
    803          ENDIF 
     974         IF ( zext(1) == obfillflt ) THEN  
     975            ! If the observation value is a fill value, set QC flag to bad  
     976            surfdataqc%nqc(jobs) = 4  
     977         ENDIF  
     978#endif defined key_cice 
    804979 
    805980      END DO 
    … …  
    814989         & zmask, & 
    815990         & zsurf, & 
    816          & zsurftmp, & 
    817          & zglamf, & 
    818          & zgphif, & 
    819          & igrdip1,& 
    820          & igrdjp1 & 
     991         & zsurftmp & 
    821992         & ) 
    822993 
    823       ! At the end of the day also deallocate night-time mean array 
    824       IF ( idayend == 0 .AND. ldnightav ) THEN 
    825          DEALLOCATE( & 
    826             & zsurfm  & 
     994      IF ( k2dint > 4 ) THEN  
     995         DEALLOCATE( &       
     996            & zglamf, &  
     997            & zgphif, &  
     998            & igrdip1,&  
     999            & igrdjp1 &  
    8271000            & ) 
    8281001      ENDIF 
    8291002 
    830       surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     1003      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )  
     1004 
     1005      ! At the end of the day also deallocate night-time mean array  
     1006      IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN 
     1007         DEALLOCATE( &  
     1008            & zsurfm  &  
     1009            & )  
     1010      ENDIF 
     1011 
     1012      IF ( kvar == surfdataqc%nvar ) THEN  
     1013         surfdataqc%nsurfup = surfdataqc%nsurfup + isurf  
     1014      ENDIF 
    8311015 
    8321016   END SUBROUTINE obs_surf_opt 

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

    r11203 r15764  
    1616   USE par_kind, ONLY : & ! Precision variables 
    1717      & wp    
     18   USE dom_oce            ! ocean space and time domain 
    1819   USE in_out_manager     ! I/O manager 
    1920   USE obs_profiles_def   ! Definitions for storage arrays for profiles 
    … …  
    5253CONTAINS 
    5354 
    54    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
    55                             kqc_cutoff ) 
     55   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, &  
     56      &                     kpi, kpj, &     
     57      &                     zmask, pglam, pgphi, &  
     58      &                     ld_nea, ld_bound_reject, &  
     59      &                     ld_seaicetypes, kqc_cutoff ) 
    5660      !!---------------------------------------------------------------------- 
    5761      !!                    ***  ROUTINE obs_pre_sla  *** 
    … …  
    8185      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data 
    8286      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening 
     87      INTEGER, INTENT(IN) :: kpi, kpj              ! Local domain sizes        
     88      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: &  
     89         & zmask        
     90      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: &  
     91         & pglam, &  
     92         & pgphi 
    8393      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
    8494      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     95      LOGICAL, INTENT(IN) :: ld_seaicetypes        ! Switch to indicate sea ice data 
    8596      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    8697      !! * Local declarations 
    … …  
    93104      INTEGER :: icycle       ! Current assimilation cycle 
    94105                              ! Counters for observations that 
    95       INTEGER :: iotdobs      !  - outside time domain 
    96       INTEGER :: iosdsobs     !  - outside space domain 
    97       INTEGER :: ilansobs     !  - within a model land cell 
    98       INTEGER :: inlasobs     !  - close to land 
    99       INTEGER :: igrdobs      !  - fail the grid search 
    100       INTEGER :: ibdysobs     !  - close to open boundary 
    101                               ! Global counters for observations that 
    102       INTEGER :: iotdobsmpp     !  - outside time domain 
    103       INTEGER :: iosdsobsmpp    !  - outside space domain 
    104       INTEGER :: ilansobsmpp    !  - within a model land cell 
    105       INTEGER :: inlasobsmpp    !  - close to land 
    106       INTEGER :: igrdobsmpp     !  - fail the grid search 
    107       INTEGER :: ibdysobsmpp  !  - close to open boundary 
    108       LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    109          & llvalid            ! SLA data selection 
    110       INTEGER :: jobs         ! Obs. loop variable 
    111       INTEGER :: jstp         ! Time loop variable 
    112       INTEGER :: inrc         ! Time index variable 
     106      INTEGER                           :: iotdobs      !  - outside time domain  
     107      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobs     !  - outside space domain  
     108      INTEGER, DIMENSION(surfdata%nvar) :: ilansobs     !  - within a model land cell  
     109      INTEGER, DIMENSION(surfdata%nvar) :: inlasobs     !  - close to land  
     110      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobs     !  - close to open boundary  
     111      INTEGER                           :: igrdobs      !  - fail the grid search        
     112                                                        ! Global counters for observations that  
     113      INTEGER                           :: iotdobsmpp   !  - outside time domain  
     114      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobsmpp  !  - outside space domain  
     115      INTEGER, DIMENSION(surfdata%nvar) :: ilansobsmpp  !  - within a model land cell  
     116      INTEGER, DIMENSION(surfdata%nvar) :: inlasobsmpp  !  - close to land  
     117      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobsmpp  !  - close to open boundary  
     118      INTEGER                           :: igrdobsmpp   !  - fail the grid search  
     119 
     120      LOGICAL, DIMENSION(:), ALLOCATABLE :: &   
     121         & llvalid            ! SLA data selection  
     122      INTEGER :: jobs         ! Obs. loop counter  
     123      INTEGER :: jvar         ! Variable loop counter  
     124      INTEGER :: jstp         ! Time loop variable  
     125      INTEGER :: inrc         ! Time index variable  
     126      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line  
    113127 
    114128      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 
    … …  
    140154      ! ----------------------------------------------------------------------- 
    141155 
    142       CALL obs_coo_tim( icycle, & 
    143          &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    144          &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
    145          &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
    146          &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
    147  
    148       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
    149        
    150       ! ----------------------------------------------------------------------- 
    151       ! Check for surface data failing the grid search 
    152       ! ----------------------------------------------------------------------- 
    153  
    154       CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, & 
    155          &              surfdata%nqc,     igrdobs                         ) 
     156      IF ( ld_seaicetypes ) THEN  
     157         CALL obs_coo_tim( icycle, &  
     158            &              iyea0,   imon0,   iday0,   ihou0, imin0,      &  
     159            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &  
     160            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &  
     161            &              surfdata%nqc,     surfdata%mstp, iotdobs,       &  
     162            &              ld_seaicetypes = ld_seaicetypes )  
     163      ELSE  
     164         CALL obs_coo_tim( icycle, &  
     165            &              iyea0,   imon0,   iday0,   ihou0, imin0,      &  
     166            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &  
     167            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &  
     168            &              surfdata%nqc,     surfdata%mstp, iotdobs        )  
     169      ENDIF  
     170 
     171      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )  
     172 
     173      ! -----------------------------------------------------------------------  
     174      ! Check for surface data failing the grid search  
     175      ! -----------------------------------------------------------------------  
     176      DO jvar = 1, surfdata%nvar  
     177         CALL obs_coo_grd( surfdata%nsurf, surfdata%mi(:,jvar), surfdata%mj(:,jvar), &  
     178            &              surfdata%nqc,     igrdobs )  
     179      END DO  
    156180 
    157181      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    … …  
    161185      ! ----------------------------------------------------------------------- 
    162186 
    163       CALL obs_coo_spc_2d( surfdata%nsurf,              & 
    164          &                 jpi,          jpj,          & 
    165          &                 surfdata%mi,   surfdata%mj,   &  
    166          &                 surfdata%rlam, surfdata%rphi, & 
    167          &                 glamt,        gphit,        & 
    168          &                 tmask(:,:,1), surfdata%nqc,  & 
    169          &                 iosdsobs,     ilansobs,     & 
    170          &                 inlasobs,     ld_nea,       & 
    171          &                 ibdysobs,     ld_bound_reject, & 
    172          &                 iqc_cutoff                     ) 
    173  
    174       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    175       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    176       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
    177       CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
     187      DO jvar = 1, surfdata%nvar  
     188         CALL obs_coo_spc_2d( surfdata%nsurf,                     &  
     189            &                 jpi,          jpj,                  &  
     190            &                 surfdata%mi(:,jvar), surfdata%mj(:,jvar), &   
     191            &                 surfdata%rlam, surfdata%rphi,       &  
     192            &                 pglam(:,:,jvar), pgphi(:,:,jvar),   &  
     193            &                 zmask(:,:,jvar), surfdata%nqc(:),   &  
     194            &                 iosdsobs(jvar), ilansobs(jvar),     &  
     195            &                 inlasobs(jvar),     ld_nea,         &  
     196            &                 ibdysobs(jvar), ld_bound_reject,    &  
     197            &                 iqc_cutoff                     )  
     198         CALL obs_mpp_sum_integer( iosdsobs(jvar), iosdsobsmpp(jvar) )  
     199         CALL obs_mpp_sum_integer( ilansobs(jvar), ilansobsmpp(jvar) )  
     200         CALL obs_mpp_sum_integer( inlasobs(jvar), inlasobsmpp(jvar) )  
     201         CALL obs_mpp_sum_integer( ibdysobs(jvar), ibdysobsmpp(jvar) )  
     202      END DO  
    178203 
    179204      ! ----------------------------------------------------------------------- 
    … …  
    204229      ! Update the total observation counter array 
    205230       
    206       IF(lwp) THEN 
    207          WRITE(numout,*) 
    208          WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', & 
    209             &            iotdobsmpp 
    210          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', & 
    211             &            igrdobsmpp 
    212          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain       = ', & 
    213             &            iosdsobsmpp 
    214          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points             = ', & 
    215             &            ilansobsmpp 
    216          IF (ld_nea) THEN 
    217             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', & 
    218                &            inlasobsmpp 
    219          ELSE 
    220             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept)    = ', & 
    221                &            inlasobsmpp 
    222          ENDIF 
    223          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
    224             &            ibdysobsmpp   
    225          WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    226             &            surfdataqc%nsurfmpp 
    227  
    228          WRITE(numout,*) 
    229          WRITE(numout,*) ' Number of observations per time step :' 
    230          WRITE(numout,*) 
    231          WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1) 
    232          WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 
    233          CALL FLUSH(numout) 
     231      IF(lwp) THEN  
     232 
     233         DO jvar = 1, surfdataqc%nvar         
     234            IF ( jvar == 1 ) THEN  
     235               cout1=TRIM(surfdataqc%cvars(1))                    
     236            ELSE  
     237               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdataqc%cvars(jvar))              
     238            ENDIF  
     239         END DO  
     240 
     241         WRITE(numout,*)  
     242         WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain                  = ', &  
     243            &            iotdobsmpp  
     244         WRITE(numout,*) ' Remaining '//TRIM(cout1)//' data that failed grid search    = ', &  
     245            &            igrdobsmpp  
     246 
     247         DO jvar = 1, surfdataqc%nvar              
     248            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data outside space domain       = ', &  
     249                &            iosdsobsmpp(jvar)  
     250             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data at land points             = ', &  
     251                &            ilansobsmpp(jvar)  
     252             IF (ld_nea) THEN  
     253                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (removed) = ', &  
     254                   &            inlasobsmpp(jvar)  
     255             ELSE  
     256                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (kept)    = ', &  
     257                   &            inlasobsmpp(jvar)  
     258             ENDIF       
     259             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near open boundary (removed) = ', &  
     260                &            ibdysobsmpp(jvar)  
     261          END DO  
     262          WRITE(numout,*) ' '//TRIM(cout1)//' data accepted = ', &  
     263             &            surfdataqc%nsurfmpp  
     264 
     265         WRITE(numout,*)  
     266         WRITE(numout,*) ' Number of observations per time step:'  
     267         WRITE(numout,*)  
     268         WRITE(numout,'(10X,A,10X,A)')'Time step',TRIM(cout1)  
     269         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'  
     270         CALL FLUSH(numout)  
    234271      ENDIF 
    235272       
    … …  
    434471 
    435472      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    436          CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
     473         CALL obs_uv_rej_pro( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    437474         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    438475         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
    … …  
    558595      &                    kobsno,                                        & 
    559596      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    560       &                    kobsqc,  kobsstp, kotdobs                      ) 
     597      &                    kobsqc,  kobsstp, kotdobs, ld_seaicetypes      ) 
    561598      !!---------------------------------------------------------------------- 
    562599      !!                    ***  ROUTINE obs_coo_tim *** 
    … …  
    606643         & kobsstp          ! Number of time steps up to the  
    607644                            ! observation time 
     645      LOGICAL, OPTIONAL, INTENT(IN) :: ld_seaicetypes 
    608646 
    609647      !! * Local declarations 
    … …  
    620658      INTEGER :: iskip 
    621659      INTEGER :: idaystp 
     660      INTEGER :: icecount 
    622661      REAL(KIND=wp) :: zminstp 
    623662      REAL(KIND=wp) :: zhoustp 
    … …  
    713752            kotdobs = kotdobs + 1 
    714753            CYCLE 
     754         ENDIF 
     755 
     756       ! Flag sea ice observations falling on initial timestep  
     757         IF ( PRESENT(ld_seaicetypes) ) THEN  
     758 
     759              IF ( ( kobsstp(jobs) == (nit000 - 1) ) ) THEN  
     760                 IF (lwp) WRITE(numout,*)( 'Sea-ice not initialised on zeroth '// &  
     761                           &    'time-step but SIT observation valid then, flagging '// &  
     762                                'in time check subroutine obs_coo_tim.' )  
     763                 kobsqc(jobs) = IBSET(kobsqc(jobs),13)  
     764                 kotdobs      = kotdobs + 1  
     765                 CYCLE  
     766              ENDIF  
    715767         ENDIF 
    716768 
    … …  
    11591211#endif  
    11601212      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
     1213         & zgdept, & 
    11611214         & zgdepw          
    11621215      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
    11631216         & zglam, &           ! Model longitude at grid points 
    1164          & zgphi              ! Model latitude at grid points 
     1217         & zgphi, &           ! Model latitude at grid points  
     1218         & zbathy             ! Index of deepest wet level at grid points 
    11651219      INTEGER, DIMENSION(2,2,kprofno) :: & 
    11661220         & igrdi, &           ! Grid i,j 
    … …  
    11701224      INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    11711225      INTEGER :: jobs, jobsp, jk, ji, jj 
     1226      REAL(KIND=wp) :: maxdepw 
    11721227 
    11731228      ! Get grid point indices 
    … …  
    12201275      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    12211276      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
     1277      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbathy), zbathy ) 
    12221278      ! Need to know the bathy depth for each observation for sco 
    12231279      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), & 
    12241280         &                  zgdepw ) 
     1281      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdept(:,:,:), &  
     1282         &                  zgdept ) 
    12251283 
    12261284      DO jobs = 1, kprofno 
    … …  
    12581316         DO jobsp = kpstart(jobs), kpend(jobs) 
    12591317 
     1318            ! Calculate max T and W depths of 2x2 grid  
     1319            maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs)  
     1320            DO jj = 1, 2  
     1321               DO ji = 1, 2  
     1322                  IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN  
     1323                     maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs)  
     1324                  END IF  
     1325               END DO  
     1326            END DO 
     1327 
    12601328            ! Flag if the observation falls outside the model spatial domain 
    12611329            IF (       ( pobslam(jobs) < -180.         )       & 
    … …  
    12641332               &  .OR. ( pobsphi(jobs) >   90.         )       & 
    12651333               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    1266                &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
     1334               &  .OR. ( pobsdep(jobsp) >= maxdepw ) ) THEN 
    12671335               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    12681336               kosdobs = kosdobs + 1 
    … …  
    13271395            ENDIF 
    13281396             
    1329             ! Set observation depth equal to that of the first model depth 
    1330             IF ( pobsdep(jobsp) <= pdep(1) ) THEN 
    1331                pobsdep(jobsp) = pdep(1)   
    1332             ENDIF 
    1333              
    13341397#if defined key_bdy 
    13351398            ! Flag if the observation falls close to the boundary rim 
    … …  
    14051468   END SUBROUTINE obs_pro_rej 
    14061469 
    1407    SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
    1408       !!---------------------------------------------------------------------- 
    1409       !!                    ***  ROUTINE obs_uv_rej *** 
     1470   SUBROUTINE obs_uv_rej_pro( profdata, knumu, knumv, kqc_cutoff ) 
     1471      !!---------------------------------------------------------------------- 
     1472      !!                    ***  ROUTINE obs_uv_rej_pro *** 
    14101473      !! 
    14111474      !! ** Purpose : Reject u if v is rejected and vice versa 
    … …  
    14391502            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN 
    14401503 
    1441             CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej') 
     1504            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej_pro') 
    14421505            RETURN 
    14431506 
    … …  
    14611524      END DO 
    14621525 
    1463    END SUBROUTINE obs_uv_rej 
     1526   END SUBROUTINE obs_uv_rej_pro 
    14641527 
    14651528END MODULE obs_prep 

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

    r11203 r15764  
    7272         & vdep,  &       !: Depth coordinate of profile data 
    7373         & vobs,  &       !: Profile data 
    74          & vmod           !: Model counterpart of the profile data vector 
     74         & vmod,  &       !: Model counterpart of the profile data vector  
     75         & vclm           !: Climatological counterpart of the profile data vector 
    7576 
    7677      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    … …  
    101102      INTEGER :: npk 
    102103      INTEGER :: nprofup  !: Observation counter used in obs_oper 
     104 
     105      LOGICAL :: lclim    !: Climatology will be calculated for this structure 
    103106 
    104107      ! Bookkeeping arrays with sizes equal to number of variables 
    … …  
    198201    
    199202   SUBROUTINE obs_prof_alloc( prof,  kvar, kext, kprof,  & 
    200       &                       ko3dt, kstp, kpi, kpj, kpk ) 
     203      &                       ko3dt, kstp, kpi, kpj, kpk, ldclim ) 
    201204      !!---------------------------------------------------------------------- 
    202205      !!                     ***  ROUTINE obs_prof_alloc  *** 
    … …  
    221224      INTEGER, INTENT(IN) :: kpj 
    222225      INTEGER, INTENT(IN) :: kpk 
     226      LOGICAL, INTENT(IN) :: ldclim 
    223227 
    224228      !!* Local variables 
    … …  
    236240      prof%npj       = kpj 
    237241      prof%npk       = kpk 
     242 
     243      prof%lclim     = ldclim 
    238244 
    239245      ! Allocate arrays of size number of variables 
    … …  
    503509            & ) 
    504510      ENDIF 
     511      IF (prof%lclim) THEN  
     512         ALLOCATE( &   
     513            & prof%var(kvar)%vclm(kobs) &  
     514            & )  
     515      ENDIF 
    505516 
    506517   END SUBROUTINE obs_prof_alloc_var 
    … …  
    537548         DEALLOCATE( &  
    538549            & prof%var(kvar)%vext  & 
     550            & )  
     551      ENDIF  
     552      IF (prof%lclim) THEN  
     553         DEALLOCATE( &   
     554            & prof%var(kvar)%vclm  & 
    539555            & ) 
    540556      ENDIF 
    … …  
    630646            &                 inprof,    invpro,    & 
    631647            &                 prof%nstp, prof%npi,  & 
    632             &                 prof%npj,  prof%npk ) 
     648            &                 prof%npj,  prof%npk,  &  
     649            &                 prof%lclim ) 
    633650      ENDIF 
    634651 
    … …  
    745762                           &                      prof%var(jvar)%vext(jj,jext) 
    746763                     END DO 
     764                     IF (newprof%lclim) THEN  
     765                        newprof%var(jvar)%vclm(invpro(jvar))   = &  
     766                           &                           prof%var(jvar)%vclm(jj)  
     767                     ENDIF 
    747768                   
    748769                     ! nvind is the index of the original variable data 
    … …  
    870891                     &                        prof%var(jvar)%vext(jj,jext) 
    871892               END DO 
     893               IF (prof%lclim) THEN  
     894                  oldprof%var(jvar)%vclm(jl)   = prof%var(jvar)%vclm(jj)  
     895               ENDIF 
    872896                
    873897            END DO 

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

    r11203 r15764  
    4646      &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    4747      &                     ldvar, ldignmis, ldsatt, & 
    48       &                     ldmod, kdailyavtypes ) 
     48      &                     ldmod, ldclim, cdvars, kdailyavtypes ) 
    4949      !!--------------------------------------------------------------------- 
    5050      !! 
    … …  
    7878      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points 
    7979      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data 
     80      LOGICAL, INTENT(IN) :: ldclim     ! Set flag to show climatology will be output 
    8081      REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS 
    8182      REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS 
     83      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    8284      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    8385         & kdailyavtypes                ! Types of daily average observations 
    … …  
    8688      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    8789      CHARACTER(len=8) :: clrefdate 
    88       CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
     90      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 
    8991      INTEGER :: jvar 
    9092      INTEGER :: ji 
    … …  
    220222            ENDIF 
    221223 
    222             IF ( jj == 1 ) THEN 
    223                ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 
    224                DO ji = 1, inpfiles(jj)%nvar 
    225                  clvars(ji) = inpfiles(jj)%cname(ji) 
    226                END DO 
     224            IF ( jj == 1 ) THEN  
     225               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )  
     226               DO ji = 1, inpfiles(jj)%nvar  
     227                 clvarsin(ji) = inpfiles(jj)%cname(ji)  
     228                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN  
     229                    CALL ctl_stop( 'Feedback file variables do not match', &  
     230                        &           ' expected variable names for this type' )  
     231                 ENDIF  
     232               END DO  
    227233            ELSE 
    228                DO ji = 1, inpfiles(jj)%nvar 
    229                   IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN 
    230                      CALL ctl_stop( 'Feedback file variables not consistent', & 
    231                         &           ' with previous files for this type' ) 
    232                   ENDIF 
    233                END DO 
     234               DO ji = 1, inpfiles(jj)%nvar  
     235                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN  
     236                     CALL ctl_stop( 'Feedback file variables not consistent', &  
     237                        &           ' with previous files for this type' )  
     238                  ENDIF  
     239               END DO  
    234240            ENDIF 
    235241 
    … …  
    500506      ENDIF 
    501507      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 
    502          &                 kstp, jpi, jpj, jpk ) 
     508         &                 kstp, jpi, jpj, jpk, ldclim ) 
    503509 
    504510      ! * Read obs/positions, QC, all variable and assign to profdata 
    … …  
    506512      profdata%nprof     = 0 
    507513      profdata%nvprot(:) = 0 
    508       profdata%cvars(:)  = clvars(:) 
     514      profdata%cvars(:)  = clvarsin(:) 
    509515      iprof = 0 
    510516 
    … …  
    696702                              profdata%var(jvar)%vmod(ivart(jvar)) = & 
    697703                                 &                inpfiles(jj)%padd(ij,ji,1,jvar) 
     704                           ENDIF 
     705                           IF ( profdata%lclim ) THEN  
     706                              profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi  
    698707                           ENDIF 
    699708                           ! Count number of profile var1 data as function of type 
    … …  
    805814      ! Deallocate temporary data 
    806815      !----------------------------------------------------------------------- 
    807       DEALLOCATE( ifileidx, iprofidx, zdat, clvars ) 
     816      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin ) 
    808817 
    809818      !----------------------------------------------------------------------- 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r11202 r15764  
    3939 
    4040   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, & 
    41       &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    42       &                     ldignmis, ldmod, ldnightav ) 
     41      &                     kvars, kextr, kstp, ddobsini, ddobsend, MeanPeriodHours, & 
     42      &                     ldignmis, ldmod, ldnightav, ldclim, ln_time_mean_sla_bkg, cdvars ) 
    4343      !!--------------------------------------------------------------------- 
    4444      !! 
    … …  
    7171      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data 
    7272      LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average 
     73      LOGICAL, INTENT(IN) :: ldclim     ! Will include climatology at obs points. 
     74      LOGICAL, INTENT(IN) :: ln_time_mean_sla_bkg     ! Will reset times to end of averaging period. 
    7375      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS 
    7476      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS 
     77      REAL(wp), INTENT(IN) :: MeanPeriodHours ! Averaging period in hours 
     78      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    7579 
    7680      !! * Local declarations 
    7781      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 
    7882      CHARACTER(len=8) :: clrefdate 
    79       CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
     83      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 
    8084      INTEGER :: ji 
    8185      INTEGER :: jj 
    8286      INTEGER :: jk 
     87      INTEGER :: jvar 
    8388      INTEGER :: iflag 
    8489      INTEGER :: inobf 
    … …  
    103108         & ityp, & 
    104109         & itypmpp 
    105       INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     110      INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 
    106111         & iobsi,    & 
    107112         & iobsj,    & 
    108          & iproc,    & 
     113         & iproc 
     114      INTEGER, DIMENSION(:), ALLOCATABLE :: &          
    109115         & iindx,    & 
    110116         & ifileidx, & 
    … …  
    122128      TYPE(obfbdata), POINTER, DIMENSION(:) :: & 
    123129         & inpfiles 
    124  
     130      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
     131       
    125132      ! Local initialization 
    126133      iobs  = 0 
    … …  
    181188               &                ldgrid = .TRUE. ) 
    182189 
     190            IF ( inpfiles(jj)%nvar /= kvars ) THEN 
     191               CALL ctl_stop( 'Feedback format error: ', & 
     192                  &           ' unexpected number of vars in feedback file' ) 
     193            ENDIF 
     194 
    183195            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 
    184196               CALL ctl_stop( 'Model not in input data' ) 
    … …  
    187199 
    188200            IF ( jj == 1 ) THEN 
    189                ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 
     201               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 
    190202               DO ji = 1, inpfiles(jj)%nvar 
    191                  clvars(ji) = inpfiles(jj)%cname(ji) 
     203                 clvarsin(ji) = inpfiles(jj)%cname(ji) 
     204                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 
     205                    CALL ctl_stop( 'Feedback file variables do not match', & 
     206                        &           ' expected variable names for this type' ) 
     207                 ENDIF 
    192208               END DO 
    193209            ELSE 
    194210               DO ji = 1, inpfiles(jj)%nvar 
    195                   IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN 
     211                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 
    196212                     CALL ctl_stop( 'Feedback file variables not consistent', & 
    197213                        &           ' with previous files for this type' ) 
    … …  
    265281 
    266282            IF ( inpfiles(jj)%nobs > 0 ) THEN 
    267                inpfiles(jj)%iproc = -1 
    268                inpfiles(jj)%iobsi = -1 
    269                inpfiles(jj)%iobsj = -1 
    270             ENDIF 
     283               inpfiles(jj)%iproc(:,:) = -1 
     284               inpfiles(jj)%iobsi(:,:) = -1 
     285               inpfiles(jj)%iobsj(:,:) = -1 
     286            ENDIF 
     287 
     288            !If SLA observations are representing a time mean then set the time 
     289            !of the obs to the end of that meaning period relative to the start of the run 
     290            IF ( ln_time_mean_sla_bkg .AND. ( TRIM( clvarsin(1) ) == 'SLA' ) ) THEN 
     291               DO ji = 1, inpfiles(jj)%nobs 
     292                  ! Only do this for obs within time window 
     293                  IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 
     294                     & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 
     295                     inpfiles(jj)%ptim(ji) = & 
     296                           & djulini(jj) + (MeanPeriodHours/24.) 
     297                  ENDIF       
     298               END DO 
     299            ENDIF    
     300             
    271301            inowin = 0 
    272302            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    278308            ALLOCATE( zlam(inowin)  ) 
    279309            ALLOCATE( zphi(inowin)  ) 
    280             ALLOCATE( iobsi(inowin) ) 
    281             ALLOCATE( iobsj(inowin) ) 
    282             ALLOCATE( iproc(inowin) ) 
     310            ALLOCATE( iobsi(inowin,kvars) ) 
     311            ALLOCATE( iobsj(inowin,kvars) ) 
     312            ALLOCATE( iproc(inowin,kvars) ) 
    283313            inowin = 0 
    284314            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    291321            END DO 
    292322 
    293             CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' ) 
     323            ! Assume anything other than velocity is on T grid 
     324            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
     325               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     326                  &                  iproc(:,1), 'U' ) 
     327               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 
     328                  &                  iproc(:,2), 'V' ) 
     329            ELSE 
     330               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     331                  &                  iproc(:,1), 'T' ) 
     332               IF ( kvars > 1 ) THEN 
     333                  DO jvar = 2, kvars 
     334                     iobsi(:,jvar) = iobsi(:,1) 
     335                     iobsj(:,jvar) = iobsj(:,1) 
     336                     iproc(:,jvar) = iproc(:,1) 
     337                  END DO 
     338               ENDIF 
     339            ENDIF 
    294340 
    295341            inowin = 0 
    … …  
    298344                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    299345                  inowin = inowin + 1 
    300                   inpfiles(jj)%iproc(ji,1) = iproc(inowin) 
    301                   inpfiles(jj)%iobsi(ji,1) = iobsi(inowin) 
    302                   inpfiles(jj)%iobsj(ji,1) = iobsj(inowin) 
     346                  DO jvar = 1, kvars 
     347                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 
     348                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 
     349                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 
     350                  END DO 
    303351               ENDIF 
    304352            END DO 
    … …  
    359407         &               iindx   ) 
    360408 
    361       CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj ) 
     409      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj, ldclim ) 
    362410 
    363411      ! Read obs/positions, QC, all variable and assign to surfdata 
    364412 
    365413      iobs = 0 
    366       surfdata%cvars(:)  = clvars(:) 
     414      surfdata%cvars(:)  = clvarsin(:) 
    367415      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'SLA' ) ) THEN 
    368416         surfdata%cext(1) = 'SSH' 
    369417         surfdata%cext(2) = 'MDT' 
     418      ENDIF 
     419      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'FBD' ) ) THEN 
     420           surfdata%cext(1) = 'freeboard' 
     421           surfdata%cext(2) = 'thick_s' 
    370422      ENDIF 
    371423      IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD' 
    … …  
    417469 
    418470               ! Coordinate search parameters 
    419                surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1) 
    420                surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1) 
    421  
     471               DO jvar = 1, kvars 
     472                  surfdata%mi(iobs,jvar) = inpfiles(jj)%iobsi(ji,jvar) 
     473                  surfdata%mj(iobs,jvar) = inpfiles(jj)%iobsj(ji,jvar) 
     474               END DO 
     475                
    422476               ! WMO number 
    423477               surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji) 
    … …  
    449503 
    450504               ! Observed value 
    451                surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
     505               DO jvar = 1, kvars                
     506                  surfdata%robs(iobs,jvar) = inpfiles(jj)%pob(1,ji,jvar) 
     507               END DO 
     508               IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN 
     509                   surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
     510                   surfdata%rext(iobs,2) = fbrmdi 
     511               ENDIF 
    452512 
    453513               ! Model and MDT is set to fbrmdi unless read from file 
    454514               IF ( ldmod ) THEN 
    455                   surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1) 
     515                  DO jvar = 1, kvars                               
     516                     surfdata%rmod(iobs,jvar) = inpfiles(jj)%padd(1,ji,1,jvar) 
     517                  END DO 
    456518                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN 
    457519                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1) 
    … …  
    459521                  ENDIF 
    460522                ELSE 
    461                   surfdata%rmod(iobs,1) = fbrmdi 
     523                  DO jvar = 1, kvars                 
     524                     surfdata%rmod(iobs,jvar) = fbrmdi 
     525                  END DO 
    462526                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi 
    463527               ENDIF 
    464  
     528                
     529               ! Initialise climatology if set 
     530               IF ( surfdata%lclim ) THEN 
     531                  DO jvar = 1, kvars 
     532                     surfdata%rclm(iobs,jvar) = fbrmdi 
     533                  END DO 
     534               ENDIF 
     535                
    465536               ! STD (obs error standard deviation) read from file and passed through obs operator 
    466537               IF ( iadd_std(jj) /= -1 ) THEN 
    … …  
    483554      !----------------------------------------------------------------------- 
    484555      IF (lwp) THEN 
    485  
     556         DO jvar = 1, surfdata%nvar        
     557            IF ( jvar == 1 ) THEN 
     558               cout1=TRIM(surfdata%cvars(1))                   
     559            ELSE 
     560               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdata%cvars(jvar))             
     561            ENDIF 
     562         END DO 
     563  
    486564         WRITE(numout,*) 
    487          WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data' 
     565         WRITE(numout,'(1X,A)')TRIM( cout1 )//' data' 
    488566         WRITE(numout,'(1X,A)')'--------------' 
    489567         DO jj = 1,8 
    … …  
    495573            & '---------------------------------------------------------------' 
    496574         WRITE(numout,'(1X,A,I8)') & 
    497             & 'Total data for variable '//TRIM( surfdata%cvars(1) )// & 
     575            & 'Total data for variable '//TRIM( cout1 )// & 
    498576            & '           = ', iobsmpp 
    499577         WRITE(numout,'(1X,A)') & 
    … …  
    506584      ! Deallocate temporary data 
    507585      !----------------------------------------------------------------------- 
    508       DEALLOCATE( ifileidx, isurfidx, zdat, clvars ) 
     586      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin ) 
    509587 
    510588      !----------------------------------------------------------------------- 

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

    r11203 r15764  
    66 
    77   !!---------------------------------------------------------------------- 
    8    !!   obs_rotvel : Rotate velocity data into N-S,E-W directorions 
     8   !!   obs_rotvel_pro :  Rotate profile velocity data into N-S,E-W directions  
     9   !!   obs_rotvel_surf : Rotate surface velocity data into N-S,E-W directions 
    910   !!---------------------------------------------------------------------- 
    1011   !! * Modules used    
    … …  
    1718   USE obs_utils                ! For error handling 
    1819   USE obs_profiles_def         ! Profile definitions 
     20   USE obs_surf_def             ! Surface definitions 
    1921   USE obs_inter_h2d            ! Horizontal interpolation 
    2022   USE obs_inter_sup            ! MPP support routines for interpolation 
    … …  
    2729   PRIVATE 
    2830 
    29    PUBLIC obs_rotvel            ! Rotate the observations 
     31   PUBLIC obs_rotvel_pro, &     ! Rotate the profile velocity observations  
     32      &   obs_rotvel_surf       ! Rotate the surface velocity observations 
    3033 
    3134   !!---------------------------------------------------------------------- 
    … …  
    3740CONTAINS 
    3841 
    39    SUBROUTINE obs_rotvel( profdata, k2dint, pu, pv ) 
     42   SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv ) 
    4043      !!--------------------------------------------------------------------- 
    4144      !! 
    … …  
    228231      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
    229232 
    230    END SUBROUTINE obs_rotvel 
     233   END SUBROUTINE obs_rotvel_pro 
     234 
     235   SUBROUTINE obs_rotvel_surf( surfdata, k2dint, pu, pv )  
     236      !!---------------------------------------------------------------------  
     237      !!  
     238      !!                   *** ROUTINE obs_rotvel_surf ***  
     239      !!  
     240      !! ** Purpose : Rotate surface velocity data into N-S,E-W directorions  
     241      !!  
     242      !! ** Method  : Interpolation of geo2ocean coefficients on U,V grid  
     243      !!              to observation point followed by a similar computations  
     244      !!              as in geo2ocean.  
     245      !!  
     246      !! ** Action  : Review if there is a better way to do this.  
     247      !!  
     248      !! References :   
     249      !!  
     250      !! History :    
     251      !!      ! :  2009-02 (K. Mogensen) : New routine  
     252      !!----------------------------------------------------------------------  
     253      !! * Modules used  
     254      !! * Arguments  
     255      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Surface data to be read  
     256      INTEGER, INTENT(IN) :: k2dint     ! Horizontal interpolation methed  
     257      REAL(wp), DIMENSION(*) :: &  
     258         & pu, &  
     259         & pv  
     260      !! * Local declarations  
     261      REAL(wp), DIMENSION(2,2,1) :: zweig  
     262      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &  
     263         & zmasku, &  
     264         & zmaskv, &  
     265         & zcoslu, &  
     266         & zsinlu, &  
     267         & zcoslv, &  
     268         & zsinlv, &  
     269         & zglamu, &  
     270         & zgphiu, &  
     271         & zglamv, &  
     272         & zgphiv  
     273      REAL(wp), DIMENSION(1) :: &  
     274         & zsinu, &  
     275         & zcosu, &  
     276         & zsinv, &  
     277         & zcosv  
     278      REAL(wp) :: zsin  
     279      REAL(wp) :: zcos  
     280      REAL(wp), DIMENSION(1) :: zobsmask  
     281      REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv  
     282      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &  
     283         & igrdiu, &  
     284         & igrdju, &  
     285         & igrdiv, &  
     286         & igrdjv  
     287      INTEGER :: ji  
     288      INTEGER :: jk  
     289 
     290      CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)   
     291 
     292      !-----------------------------------------------------------------------  
     293      ! Allocate data for message parsing and interpolation  
     294      !-----------------------------------------------------------------------  
     295 
     296      ALLOCATE( &  
     297         & igrdiu(2,2,surfdata%nsurf), &  
     298         & igrdju(2,2,surfdata%nsurf), &  
     299         & zglamu(2,2,surfdata%nsurf), &  
     300         & zgphiu(2,2,surfdata%nsurf), &  
     301         & zmasku(2,2,surfdata%nsurf), &  
     302         & zcoslu(2,2,surfdata%nsurf), &  
     303         & zsinlu(2,2,surfdata%nsurf), &  
     304         & igrdiv(2,2,surfdata%nsurf), &  
     305         & igrdjv(2,2,surfdata%nsurf), &  
     306         & zglamv(2,2,surfdata%nsurf), &  
     307         & zgphiv(2,2,surfdata%nsurf), &  
     308         & zmaskv(2,2,surfdata%nsurf), &  
     309         & zcoslv(2,2,surfdata%nsurf), &  
     310         & zsinlv(2,2,surfdata%nsurf)  &  
     311         & )  
     312 
     313      !-----------------------------------------------------------------------  
     314      ! Receive the angles on the U and V grids.  
     315      !-----------------------------------------------------------------------  
     316 
     317      CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv )  
     318 
     319      DO ji = 1, surfdata%nsurf  
     320         igrdiu(1,1,ji) = surfdata%mi(ji,1)-1  
     321         igrdju(1,1,ji) = surfdata%mj(ji,1)-1  
     322         igrdiu(1,2,ji) = surfdata%mi(ji,1)-1  
     323         igrdju(1,2,ji) = surfdata%mj(ji,1)  
     324         igrdiu(2,1,ji) = surfdata%mi(ji,1)  
     325         igrdju(2,1,ji) = surfdata%mj(ji,1)-1  
     326         igrdiu(2,2,ji) = surfdata%mi(ji,1)  
     327         igrdju(2,2,ji) = surfdata%mj(ji,1)  
     328         igrdiv(1,1,ji) = surfdata%mi(ji,2)-1  
     329         igrdjv(1,1,ji) = surfdata%mj(ji,2)-1  
     330         igrdiv(1,2,ji) = surfdata%mi(ji,2)-1  
     331         igrdjv(1,2,ji) = surfdata%mj(ji,2)  
     332         igrdiv(2,1,ji) = surfdata%mi(ji,2)  
     333         igrdjv(2,1,ji) = surfdata%mj(ji,2)-1  
     334         igrdiv(2,2,ji) = surfdata%mi(ji,2)  
     335         igrdjv(2,2,ji) = surfdata%mj(ji,2)  
     336      END DO  
     337 
     338      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &  
     339         &                  glamu, zglamu )  
     340      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &  
     341         &                  gphiu, zgphiu )  
     342      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &  
     343         &                  umask(:,:,1), zmasku )  
     344      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &  
     345         &                  zsingu, zsinlu )  
     346      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &  
     347         &                  zcosgu, zcoslu )  
     348      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &  
     349         &                  glamv, zglamv )  
     350      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &  
     351         &                  gphiv, zgphiv )  
     352      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &  
     353         &                  vmask(:,:,1), zmaskv )  
     354      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &  
     355         &                  zsingv, zsinlv )  
     356      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &  
     357         &                  zcosgv, zcoslv )  
     358 
     359      DO ji = 1, surfdata%nsurf  
     360 
     361         CALL obs_int_h2d_init( 1, 1, k2dint, &  
     362            &                   surfdata%rlam(ji), surfdata%rphi(ji), &  
     363            &                   zglamu(:,:,ji), zgphiu(:,:,ji),       &  
     364            &                   zmasku(:,:,ji), zweig, zobsmask )  
     365 
     366         CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji),  zsinu )  
     367 
     368         CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji),  zcosu )  
     369 
     370         CALL obs_int_h2d_init( 1, 1, k2dint, &  
     371            &                   surfdata%rlam(ji), surfdata%rphi(ji), &  
     372            &                   zglamv(:,:,ji), zgphiv(:,:,ji),       &  
     373            &                   zmaskv(:,:,ji), zweig, zobsmask )  
     374 
     375         CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji),  zsinv )  
     376 
     377         CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji),  zcosv )  
     378 
     379         ! Assume that the angle at observation point is the   
     380         ! mean of u and v cosines/sines  
     381 
     382         zcos = 0.5_wp * ( zcosu(1) + zcosv(1) )  
     383         zsin = 0.5_wp * ( zsinu(1) + zsinv(1) )  
     384 
     385         IF ( ( surfdata%rmod(ji,1) /= fbrmdi ) .AND. &  
     386            & ( surfdata%rmod(ji,2) /= fbrmdi ) ) THEN  
     387            pu(ji) = surfdata%rmod(ji,1) * zcos - &  
     388               &     surfdata%rmod(ji,2) * zsin  
     389            pv(ji) = surfdata%rmod(ji,2) * zcos + &  
     390               &     surfdata%rmod(ji,1) * zsin  
     391         ELSE  
     392            pu(ji) = fbrmdi  
     393            pv(ji) = fbrmdi  
     394         ENDIF  
     395 
     396      END DO  
     397 
     398      DEALLOCATE( &  
     399         & igrdiu, &  
     400         & igrdju, &  
     401         & zglamu, &  
     402         & zgphiu, &  
     403         & zmasku, &  
     404         & zcoslu, &  
     405         & zsinlu, &  
     406         & igrdiv, &  
     407         & igrdjv, &  
     408         & zglamv, &  
     409         & zgphiv, &  
     410         & zmaskv, &  
     411         & zcoslv, &  
     412         & zsinlv  &  
     413         & )  
     414 
     415      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)   
     416 
     417   END SUBROUTINE obs_rotvel_surf  
    231418 
    232419END MODULE obs_rot_vel 

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

    r11203 r15764  
    5252      INTEGER :: nrec       !: Number of surface observation records in window 
    5353 
     54      LOGICAL :: lclim      !: Climatology will be calculated for this structure 
     55 
    5456      ! Arrays with size equal to the number of surface observations 
    5557 
    5658      INTEGER, POINTER, DIMENSION(:) :: & 
    57          & mi,   &        !: i-th grid coord. for interpolating to surface observation 
    58          & mj,   &        !: j-th grid coord. for interpolating to surface observation 
    5959         & mt,   &        !: time record number for gridded data 
    6060         & nsidx,&        !: Surface observation number 
    … …  
    6969         & ntyp           !: Type of surface observation product 
    7070 
     71      INTEGER, POINTER, DIMENSION(:,:) :: &  
     72         & mi,   &        !: i-th grid coord. for interpolating to surface observation  
     73         & mj             !: j-th grid coord. for interpolating to surface observation 
     74 
    7175      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    7276         & cvars          !: Variable names 
    … …  
    8488      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    8589         & robs, &        !: Surface observation  
    86          & rmod           !: Model counterpart of the surface observation vector 
     90         & rmod, &        !: Model counterpart of the surface observation vector  
     91         & rclm           !: Climatological counterpart of the surface observation vector 
    8792 
    8893      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    8994         & rext           !: Extra fields interpolated to observation points 
    9095 
    91       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
     96      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: & 
    9297         & vdmean         !: Time averaged of model field 
    9398 
    … …  
    111116 
    112117      ! Is this a gridded product? 
    113       
    114118      LOGICAL :: lgrid 
    115119 
    … …  
    124128CONTAINS 
    125129    
    126    SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj ) 
     130   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj, ldclim ) 
    127131      !!---------------------------------------------------------------------- 
    128132      !!                     ***  ROUTINE obs_surf_alloc  *** 
    … …  
    143147      INTEGER, INTENT(IN) :: kpi     ! Number of 3D grid points 
    144148      INTEGER, INTENT(IN) :: kpj 
     149      LOGICAL, INTENT(IN) :: ldclim 
    145150 
    146151      !!* Local variables 
    … …  
    157162      surf%npi      = kpi 
    158163      surf%npj      = kpj 
     164      surf%lclim    = ldclim 
    159165 
    160166      ! Allocate arrays of size number of variables 
    … …  
    171177 
    172178      ALLOCATE( & 
    173          & surf%mi(ksurf),      & 
    174          & surf%mj(ksurf),      & 
    175179         & surf%mt(ksurf),      & 
    176180         & surf%nsidx(ksurf),   & 
    … …  
    196200 
    197201      ALLOCATE( &  
     202         & surf%mi(ksurf,kvar),   &  
     203         & surf%mj(ksurf,kvar),   & 
    198204         & surf%robs(ksurf,kvar), & 
    199205         & surf%rmod(ksurf,kvar)  & 
    200206         & )    
    201207 
     208      IF (surf%lclim) ALLOCATE( surf%rclm(ksurf,kvar) ) 
     209 
    202210      ! Allocate arrays of number of extra fields at observation points 
    203211 
    … …  
    223231 
    224232      ALLOCATE( & 
    225          & surf%vdmean(kpi,kpj) & 
     233         & surf%vdmean(kpi,kpj,kvar) & 
    226234         & ) 
    227235 
    … …  
    293301         & ) 
    294302 
     303      IF (surf%lclim) DEALLOCATE( surf%rclm ) 
    295304      ! Deallocate arrays of number of extra fields at observation points 
    296305 
    … …  
    371380      IF ( lallocate ) THEN 
    372381         CALL obs_surf_alloc( newsurf,  insurf, surf%nvar, & 
    373             & surf%nextra, surf%nstp, surf%npi, surf%npj ) 
     382            & surf%nextra, surf%nstp, surf%npi, surf%npj, surf%lclim ) 
    374383      ENDIF 
    375384 
    … …  
    418427               newsurf%robs(insurf,jk)  = surf%robs(ji,jk) 
    419428               newsurf%rmod(insurf,jk)  = surf%rmod(ji,jk) 
     429               IF (newsurf%lclim) newsurf%rclm(insurf,jk) = surf%rclm(ji,jk) 
    420430                
    421431            END DO 
    … …  
    514524            oldsurf%robs(jj,jk)  = surf%robs(ji,jk) 
    515525            oldsurf%rmod(jj,jk)  = surf%rmod(ji,jk) 
     526            IF (surf%lclim) oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk) 
    516527 
    517528         END DO 

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

    r11203 r15764  
    2828   USE obs_const 
    2929   USE obs_mpp              ! MPP support routines for observation diagnostics 
    30    USE lib_mpp        ! MPP routines 
     30   USE lib_mpp              ! MPP routines 
    3131 
    3232   IMPLICIT NONE 
    … …  
    9595      INTEGER :: je 
    9696      INTEGER :: iadd 
     97      INTEGER :: iadd_clm ! 1 if climatology present 
    9798      INTEGER :: iext 
    9899      REAL(wp) :: zpres 
     100 
     101      iadd_clm = 0   
     102      IF ( profdata%lclim ) iadd_clm = 1 
    99103 
    100104      IF ( PRESENT( padd ) ) THEN 
    … …  
    137141         fbdata%caddunit(1,1) = 'Degrees centigrade' 
    138142         fbdata%caddunit(1,2) = 'PSU' 
     143         IF ( profdata%lclim ) THEN  
     144            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature'  
     145            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity'  
     146            fbdata%caddunit(2,1) = 'Degrees centigrade'  
     147            fbdata%caddunit(2,2) = 'PSU'  
     148         ENDIF 
    139149         fbdata%cgrid(:)      = 'T' 
    140150         DO je = 1, iext 
    … …  
    170180         fbdata%caddunit(1,1) = 'm/s' 
    171181         fbdata%caddunit(1,2) = 'm/s' 
     182         IF ( profdata%lclim ) THEN  
     183            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity'  
     184            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity'  
     185            fbdata%caddunit(2,1) = 'm/s'  
     186            fbdata%caddunit(2,2) = 'm/s'  
     187         ENDIF 
    172188         fbdata%cgrid(1)      = 'U'  
    173189         fbdata%cgrid(2)      = 'V' 
    … …  
    252268         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 
    253269         fbdata%caddunit(1,1) = clunits 
     270         IF ( profdata%lclim ) THEN  
     271            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname)  
     272            fbdata%caddunit(2,1) = clunits  
     273         ENDIF 
    254274         fbdata%cgrid(:)      = clgrid 
    255275         DO je = 1, iext 
    … …  
    266286 
    267287      fbdata%caddname(1)   = 'Hx' 
     288 
     289      IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm)   = 'CLM' 
    268290 
    269291      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 
    … …  
    319341            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    320342               ik = profdata%var(jvar)%nvlidx(jk) 
    321                fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
    322343               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk) 
    323344               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk) 
    … …  
    327348                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 
    328349                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 
    329                   fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111') 
     350                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111') 
    330351               ELSE 
    331352                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) 
    … …  
    333354               ENDIF 
    334355               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk) 
    335                DO ja = 1, iadd 
    336                   fbdata%padd(ik,jo,1+ja,jvar) = & 
     356               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)  
     357               IF ( profdata%lclim ) THEN             
     358                  fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk)       
     359               ENDIF                
     360               DO ja = 1, iadd  
     361                  fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = & 
    337362                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja)) 
    338363               END DO 
    … …  
    416441      INTEGER :: ja 
    417442      INTEGER :: je 
     443      INTEGER :: jv 
    418444      INTEGER :: iadd 
    419445      INTEGER :: iext 
    420446      INTEGER :: indx_std 
    421447      INTEGER :: iadd_std 
    422  
     448      INTEGER :: iadd_clm       
     449      INTEGER :: iadd_mdt 
     450 
     451      IF ( PRESENT( pext ) ) THEN  
     452         iext = pext%inum  
     453      ELSE  
     454         iext = 0  
     455      ENDIF  
     456 
     457      ! Set up number of additional variables to be ouput:  
     458      ! Hx, CLM, STD, MDT... 
    423459      IF ( PRESENT( padd ) ) THEN 
    424460         iadd = padd%inum 
    425461      ELSE 
    426462         iadd = 0 
    427       ENDIF 
    428  
    429       IF ( PRESENT( pext ) ) THEN 
    430          iext = pext%inum 
    431       ELSE 
    432          iext = 0 
    433463      ENDIF 
    434464 
    … …  
    444474      ENDIF 
    445475       
     476      iadd_clm = 0  
     477      IF ( surfdata%lclim ) iadd_clm = 1  
     478 
     479      iadd_mdt = 0  
     480      IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1  
     481 
    446482      CALL init_obfbdata( fbdata ) 
    447483 
    … …  
    451487         ! SLA needs special treatment because of MDT, so is all done here 
    452488         ! Other variables are done more generically 
    453  
    454          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    455             &                 2 + iadd + iadd_std, 1 + iext, .TRUE. ) 
     489         ! No climatology for SLA, MDT is our best estimate of that and is already output. 
     490 
     491         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &  
     492            &                 1 + iadd_mdt + iadd_std + iadd, &  
     493            &                 1 + iext, .TRUE. ) 
    456494 
    457495         clfiletype = 'slafb' 
    … …  
    493531         clgrid     = 'T' 
    494532 
     533      CASE('SIT')  
     534 
     535         clfiletype    = 'sitfb'  
     536         cllongname    = 'Sea ice thickness'  
     537         clunits       = 'm'  
     538         clgrid        = 'T'  
     539 
     540      CASE('FBD')  
     541 
     542         clfiletype = 'fbdfb'  
     543         ! Change label from FBD ("freeboard") to SIT ("Sea Ice Thickness")  
     544         surfdata%cvars(1) = 'SIT'  
     545         cllongname = 'Sea ice thickness'  
     546         clunits    = 'm'  
     547         clgrid     = 'T'  
     548 
    495549      CASE('SSS') 
    496550 
    … …  
    500554         clgrid     = 'T' 
    501555 
    502       CASE('SLCHLTOT','LOGCHL','LogChl','logchl') 
    503  
    504          clfiletype = 'slchltotfb' 
    505          cllongname = 'Surface total log10(chlorophyll)' 
    506          clunits    = 'log10(mg/m3)' 
    507          clgrid     = 'T' 
    508  
    509       CASE('SLCHLDIA') 
    510  
    511          clfiletype = 'slchldiafb' 
    512          cllongname = 'Surface diatom log10(chlorophyll)' 
    513          clunits    = 'log10(mg/m3)' 
    514          clgrid     = 'T' 
    515  
    516       CASE('SLCHLNON') 
    517  
    518          clfiletype = 'slchlnonfb' 
    519          cllongname = 'Surface non-diatom log10(chlorophyll)' 
    520          clunits    = 'log10(mg/m3)' 
    521          clgrid     = 'T' 
    522  
    523       CASE('SLCHLDIN') 
    524  
    525          clfiletype = 'slchldinfb' 
    526          cllongname = 'Surface dinoflagellate log10(chlorophyll)' 
    527          clunits    = 'log10(mg/m3)' 
    528          clgrid     = 'T' 
    529  
    530       CASE('SLCHLMIC') 
    531  
    532          clfiletype = 'slchlmicfb' 
    533          cllongname = 'Surface microphytoplankton log10(chlorophyll)' 
    534          clunits    = 'log10(mg/m3)' 
    535          clgrid     = 'T' 
    536  
    537       CASE('SLCHLNAN') 
    538  
    539          clfiletype = 'slchlnanfb' 
    540          cllongname = 'Surface nanophytoplankton log10(chlorophyll)' 
    541          clunits    = 'log10(mg/m3)' 
    542          clgrid     = 'T' 
    543  
    544       CASE('SLCHLPIC') 
    545  
    546          clfiletype = 'slchlpicfb' 
    547          cllongname = 'Surface picophytoplankton log10(chlorophyll)' 
    548          clunits    = 'log10(mg/m3)' 
    549          clgrid     = 'T' 
    550  
    551       CASE('SCHLTOT') 
    552  
    553          clfiletype = 'schltotfb' 
    554          cllongname = 'Surface total chlorophyll' 
    555          clunits    = 'mg/m3' 
    556          clgrid     = 'T' 
    557  
    558       CASE('SLPHYTOT') 
    559  
    560          clfiletype = 'slphytotfb' 
    561          cllongname = 'Surface total log10(phytoplankton carbon)' 
    562          clunits    = 'log10(mmolC/m3)' 
    563          clgrid     = 'T' 
    564  
    565       CASE('SLPHYDIA') 
    566  
    567          clfiletype = 'slphydiafb' 
    568          cllongname = 'Surface diatom log10(phytoplankton carbon)' 
    569          clunits    = 'log10(mmolC/m3)' 
    570          clgrid     = 'T' 
    571  
    572       CASE('SLPHYNON') 
    573  
    574          clfiletype = 'slphynonfb' 
    575          cllongname = 'Surface non-diatom log10(phytoplankton carbon)' 
    576          clunits    = 'log10(mmolC/m3)' 
    577          clgrid     = 'T' 
    578  
    579       CASE('SSPM') 
    580  
    581          clfiletype = 'sspmfb' 
    582          cllongname = 'Surface suspended particulate matter' 
    583          clunits    = 'g/m3' 
    584          clgrid     = 'T' 
    585  
    586       CASE('SFCO2','FCO2','fCO2','fco2') 
    587  
    588          clfiletype = 'sfco2fb' 
    589          cllongname = 'Surface fugacity of carbon dioxide' 
    590          clunits    = 'uatm' 
    591          clgrid     = 'T' 
    592  
    593       CASE('SPCO2') 
    594  
    595          clfiletype = 'spco2fb' 
    596          cllongname = 'Surface partial pressure of carbon dioxide' 
    597          clunits    = 'uatm' 
    598          clgrid     = 'T' 
     556      CASE('UVEL')  
     557 
     558         clfiletype    = 'ssvfb'  
     559         cllongname    = 'Eastward sea surface velocity'  
     560         clunits       = 'm s-1'  
     561         clgrid        = 'U'  
     562         cllongname    = 'Northward sea surface velocity'  
     563         clunits       = 'm s-1'  
     564         clgrid        = 'V'  
     565 
     566      CASE('SLCHLTOT')  
     567 
     568         clfiletype    = 'slchltotfb'  
     569         cllongname    = 'Surface total log10(chlorophyll)'  
     570         clunits       = 'log10(mg/m3)'  
     571         clgrid        = 'T'  
     572 
     573      CASE('SLCHLDIA')  
     574 
     575         clfiletype    = 'slchldiafb'  
     576         cllongname    = 'Surface diatom log10(chlorophyll)'  
     577         clunits       = 'log10(mg/m3)'  
     578         clgrid        = 'T'  
     579 
     580      CASE('SLCHLNON')  
     581 
     582         clfiletype    = 'slchlnonfb'  
     583         cllongname    = 'Surface non-diatom log10(chlorophyll)'  
     584         clunits       = 'log10(mg/m3)'  
     585         clgrid        = 'T'  
     586 
     587      CASE('SLCHLDIN')  
     588 
     589         clfiletype    = 'slchldinfb'  
     590         cllongname    = 'Surface dinoflagellate log10(chlorophyll)'  
     591         clunits       = 'log10(mg/m3)'  
     592         clgrid        = 'T'  
     593 
     594      CASE('SLCHLMIC')  
     595 
     596         clfiletype    = 'slchlmicfb'  
     597         cllongname    = 'Surface microphytoplankton log10(chlorophyll)'  
     598         clunits       = 'log10(mg/m3)'  
     599         clgrid        = 'T'  
     600 
     601      CASE('SLCHLNAN')  
     602 
     603         clfiletype    = 'slchlnanfb'  
     604         cllongname    = 'Surface nanophytoplankton log10(chlorophyll)'  
     605         clunits       = 'log10(mg/m3)'  
     606         clgrid        = 'T'  
     607 
     608      CASE('SLCHLPIC')  
     609 
     610         clfiletype    = 'slchlpicfb'  
     611         cllongname    = 'Surface picophytoplankton log10(chlorophyll)'  
     612         clunits       = 'log10(mg/m3)'  
     613         clgrid        = 'T'  
     614 
     615      CASE('SCHLTOT')  
     616 
     617         clfiletype = 'schltotfb'  
     618         cllongname = 'Surface total chlorophyll'  
     619         clunits    = 'mg/m3'  
     620         clgrid     = 'T'  
     621 
     622      CASE('SLPHYTOT')  
     623 
     624         clfiletype    = 'slphytotfb'  
     625         cllongname    = 'Surface total log10(phytoplankton carbon)'  
     626         clunits       = 'log10(mmolC/m3)'  
     627         clgrid        = 'T'  
     628 
     629      CASE('SLPHYDIA')  
     630 
     631         clfiletype    = 'slphydiafb'  
     632         cllongname    = 'Surface diatom log10(phytoplankton carbon)'  
     633         clunits       = 'log10(mmolC/m3)'  
     634         clgrid        = 'T'  
     635 
     636      CASE('SLPHYNON')  
     637 
     638         clfiletype    = 'slphynonfb'  
     639         cllongname    = 'Surface non-diatom log10(phytoplankton carbon)'  
     640         clunits       = 'log10(mmolC/m3)'  
     641         clgrid        = 'T'  
     642 
     643      CASE('SSPM')  
     644 
     645         clfiletype    = 'sspmfb'  
     646         cllongname    = 'Surface suspended particulate matter'  
     647         clunits       = 'g/m3'  
     648         clgrid        = 'T'  
     649 
     650      CASE('SKD490')  
     651 
     652         clfiletype    = 'skd490fb'  
     653         cllongname    = 'Surface attenuation coefficient of downwelling radiation at 490 nm'  
     654         clunits       = 'm-1'  
     655         clgrid        = 'T'  
     656 
     657      CASE('SFCO2')  
     658 
     659         clfiletype    = 'sfco2fb'  
     660         cllongname    = 'Surface fugacity of carbon dioxide'  
     661         clunits       = 'uatm'  
     662         clgrid        = 'T'  
     663 
     664      CASE('SPCO2')  
     665 
     666         clfiletype    = 'spco2fb'  
     667         cllongname    = 'Surface partial pressure of carbon dioxide'  
     668         clunits       = 'uatm'  
     669         clgrid        = 'T' 
    599670 
    600671      CASE DEFAULT 
    … …  
    607678      ! Remaining variables treated more generically 
    608679 
    609       IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 
    610        
    611          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    612             &                 1 + iadd + iadd_std, iext, .TRUE. ) 
    613  
    614          fbdata%cname(1)      = surfdata%cvars(1) 
    615          fbdata%coblong(1)    = cllongname 
    616          fbdata%cobunit(1)    = clunits 
    617          DO je = 1, iext 
    618             fbdata%cextname(je) = pext%cdname(je) 
    619             fbdata%cextlong(je) = pext%cdlong(je,1) 
    620             fbdata%cextunit(je) = pext%cdunit(je,1) 
    621          END DO 
    622          IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
    623             fbdata%caddlong(1,1) = 'Model interpolated ICE' 
    624          ELSE 
    625             fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1)) 
    626          ENDIF 
    627          fbdata%caddunit(1,1) = clunits 
    628          fbdata%cgrid(1)      = clgrid 
    629          DO ja = 1, iadd 
    630             fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja) 
    631             fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1) 
    632             fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1) 
    633          END DO 
    634  
     680      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN  
     681 
     682         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, &  
     683            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )  
     684 
     685         DO jv = 1, surfdata%nvar  
     686            fbdata%cname(jv)      = surfdata%cvars(jv)  
     687            fbdata%coblong(jv)    = cllongname  
     688            fbdata%cobunit(jv)    = clunits 
     689         END DO  
     690         DO je = 1, iext  
     691            fbdata%cextname(je) = pext%cdname(je)  
     692            fbdata%cextlong(je) = pext%cdlong(je,1)  
     693            fbdata%cextunit(je) = pext%cdunit(je,1)  
     694         END DO  
     695         DO jv = 1, surfdata%nvar           
     696            IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN  
     697               fbdata%caddlong(1,jv) = 'Model interpolated ICE'  
     698            ELSE  
     699               fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv))  
     700            ENDIF  
     701            fbdata%caddunit(1,jv) = clunits 
     702            fbdata%cgrid(jv)      = clgrid 
     703         END DO              
     704         DO ja = 1, iadd  
     705            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)  
     706            DO jv = 1, surfdata%nvar                       
     707               fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv)  
     708               fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv)  
     709            END DO  
     710         END DO  
    635711      ENDIF 
    636712       
    637713      fbdata%caddname(1)   = 'Hx' 
    638       IF ( indx_std /= -1 ) THEN 
    639          IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1 
    640          fbdata%caddname(1+iadd_std)   = surfdata%cext(indx_std) 
    641          fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation' 
    642          fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1) 
     714      IF ( indx_std /= -1 ) THEN  
     715         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std)  
     716         DO jv = 1, surfdata%nvar                          
     717            fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'  
     718            fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)  
     719         END DO  
     720      ENDIF  
     721 
     722      IF ( surfdata%lclim ) THEN  
     723         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM'  
     724         DO jv = 1, surfdata%nvar                                   
     725            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology'  
     726            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1)  
     727         END DO  
    643728      ENDIF 
    644729       
    … …  
    649734         WRITE(numout,*)'obs_wri_surf :' 
    650735         WRITE(numout,*)'~~~~~~~~~~~~~' 
    651          WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname) 
     736         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname) 
    652737      ENDIF 
    653738 
    … …  
    663748            fbdata%ioqc(jo)    = 4 
    664749            fbdata%ioqcf(1,jo) = 0 
    665             fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
     750            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111') 
    666751         ELSE 
    667752            fbdata%ioqc(jo)    = surfdata%nqc(jo) 
    … …  
    675760         fbdata%kindex(jo)    = surfdata%nsfil(jo) 
    676761         IF (ln_grid_global) THEN 
    677             fbdata%iobsi(jo,1) = surfdata%mi(jo) 
    678             fbdata%iobsj(jo,1) = surfdata%mj(jo) 
     762            DO jv = 1, surfdata%nvar  
     763               fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv)  
     764               fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv)  
     765            END DO 
    679766         ELSE 
    680             fbdata%iobsi(jo,1) = mig(surfdata%mi(jo)) 
    681             fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo)) 
     767            DO jv = 1, surfdata%nvar           
     768               fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv))  
     769               fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv))  
     770            END DO 
    682771         ENDIF 
    683772         CALL greg2jul( 0, & 
    … …  
    689778            &           fbdata%ptim(jo),   & 
    690779            &           krefdate = 19500101 ) 
    691          fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1) 
    692          IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1) 
    693          fbdata%pob(1,jo,1)    = surfdata%robs(jo,1)  
    694780         fbdata%pdep(1,jo)     = 0.0 
    695781         fbdata%idqc(1,jo)     = 0 
    696782         fbdata%idqcf(:,1,jo)  = 0 
    697          IF ( surfdata%nqc(jo) > 255 ) THEN 
    698             fbdata%ivqc(jo,1)       = 4 
    699             fbdata%ivlqc(1,jo,1)    = 4 
    700             fbdata%ivlqcf(1,1,jo,1) = 0 
    701             fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
    702          ELSE 
    703             fbdata%ivqc(jo,1)       = surfdata%nqc(jo) 
    704             fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo) 
    705             fbdata%ivlqcf(:,1,jo,1) = 0 
    706          ENDIF 
    707          fbdata%iobsk(1,jo,1)  = 0 
    708          IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 
    709          IF ( indx_std /= -1 ) THEN 
    710             fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std) 
    711          ENDIF 
    712           
    713          DO ja = 1, iadd 
    714             fbdata%padd(1,jo,2+iadd_std+ja,1) = & 
    715                & surfdata%rext(jo,padd%ipoint(ja)) 
    716          END DO 
    717          DO je = 1, iext 
    718             fbdata%pext(1,jo,1+je) = & 
    719                & surfdata%rext(jo,pext%ipoint(je)) 
     783         DO jv = 1, surfdata%nvar   
     784            fbdata%pob(1,jo,jv)    = surfdata%robs(jo,jv)  
     785 
     786            IF ( surfdata%nqc(jo) > 255 ) THEN  
     787               fbdata%ivqc(jo,jv)       = 4  
     788               fbdata%ivlqc(1,jo,jv)    = 4  
     789               fbdata%ivlqcf(1,1,jo,jv) = 0  
     790               fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000000011111111')  
     791            ELSE  
     792               fbdata%ivqc(jo,jv)       = surfdata%nqc(jo)  
     793               fbdata%ivlqc(1,jo,jv)    = surfdata%nqc(jo)  
     794               fbdata%ivlqcf(:,1,jo,jv) = 0  
     795            ENDIF  
     796            fbdata%iobsk(1,jo,jv)  = 0  
     797 
     798            ! Additional variables.  
     799            ! Hx is always the first additional variable  
     800            fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv)  
     801            ! MDT is output as an additional variable if SLA obs type  
     802            IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN  
     803               fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1)  
     804            ENDIF      
     805            ! STD is output as an additional variable if available  
     806            IF ( indx_std /= -1 ) THEN  
     807               fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std)  
     808            ENDIF  
     809            ! CLM is output as an additional variable if available  
     810            IF ( surfdata%lclim ) THEN  
     811               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv)  
     812            ENDIF  
     813            ! Then other additional variables are output  
     814            DO ja = 1, iadd  
     815               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = &  
     816                  & surfdata%rext(jo,padd%ipoint(ja))  
     817            END DO  
     818         END DO           
     819         ! Extra variables  
     820         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)           
     821         DO je = 1, iext  
     822            fbdata%pext(1,jo,1+je) = &  
     823               & surfdata%rext(jo,pext%ipoint(je))  
    720824         END DO 
    721825      END DO 

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

    r11203 r15764  
    6262         z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      ) 
    6363         z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) ) 
    64          IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp 
    65  
    66          zsum = z1dm + z1dp 
    6764          
    68          IF ( k1dint == 0 ) THEN 
    69  
    70             !----------------------------------------------------------------- 
    71             !  Linear interpolation 
    72             !----------------------------------------------------------------- 
    73             pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) & 
    74                &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum 
    75  
    76          ELSEIF ( k1dint == 1 ) THEN 
    77  
    78             !----------------------------------------------------------------- 
    79             ! Cubic spline interpolation 
    80             !----------------------------------------------------------------- 
    81             zsum2 = zsum * zsum 
    82             pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) & 
    83                &           + z1dp                             * pobsk (kkco(jdep)  ) & 
    84                &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 
    85                &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) & 
    86                &             ) / 6.0_wp                                              & 
    87                &          ) / zsum 
    88  
     65         ! Where both levels are masked, return a fill value  
     66         IF ( ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .AND. (pobsmask(kkco(jdep)) == 0.0_wp) ) THEN  
     67            pobs(jdep)  = 99999.  
     68         ELSE  
     69           
     70            ! Where upper level is masked (e.g., under ice cavity), only use deeper level  
     71            ! otherwise where ob is at or above upper level model T-point,   
     72            ! use upper model level rather than extrapolate  
     73            IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) THEN  
     74               z1dm = 0.0_wp  
     75            ELSE IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) THEN  
     76               z1dp = 0.0_wp     
     77            END IF     
     78  
     79            ! Where deeper level is masked (e.g., near sea bed), only use upper level  
     80            ! otherwise where ob is at or below deeper level model T-point,   
     81            ! use deeper model level rather than extrapolate  
     82            IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN  
     83               z1dp = 0.0_wp  
     84            ELSE IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) THEN  
     85               z1dm = 0.0_wp     
     86            END IF     
     87  
     88            zsum = z1dm + z1dp  
     89           
     90            IF ( k1dint == 0 ) THEN  
     91  
     92               !-----------------------------------------------------------------  
     93               !  Linear interpolation  
     94               !-----------------------------------------------------------------  
     95               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &  
     96                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum  
     97  
     98            ELSEIF ( k1dint == 1 ) THEN  
     99  
     100               !-----------------------------------------------------------------  
     101               ! Cubic spline interpolation  
     102               !-----------------------------------------------------------------  
     103               zsum2 = zsum * zsum  
     104               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &  
     105                  &           + z1dp                             * pobsk (kkco(jdep)  ) &  
     106                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &  
     107                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &  
     108                  &             ) / 6.0_wp                                              &  
     109                  &          ) / zsum  
     110  
     111            ENDIF  
    89112         ENDIF 
    90113      END DO 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r11203 r15764  
    5656   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1) 
    5757 
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tclim    !: temperature climatology on each time step(Celcius)  
     59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sclim    !: salinity climatology on each time step (psu) 
     60 
    5861   !! * Substitutions 
    5962#  include "domzgr_substitute.h90" 
    … …  
    7073      !!                ***  FUNCTION tra_dmp_alloc  *** 
    7174      !!---------------------------------------------------------------------- 
    72       ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc ) 
     75      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), &  
     76         &      tclim(jpi,jpj,jpk) , sclim(jpi,jpj,jpk), STAT=tra_dmp_alloc ) 
    7377      ! 
    7478      IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc ) 
    … …  
    110114      !                           !==   input T-S data at kt   ==! 
    111115      CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt 
     116 
     117      tclim(:,:,:) = zts_dta(:,:,:,jp_tem)  
     118      sclim(:,:,:) = zts_dta(:,:,:,jp_sal) 
    112119      ! 
    113120      SELECT CASE ( nn_zdmp )     !==    type of damping   ==! 
Note: See TracChangeset for help on using the changeset viewer.