Changeset 12610


Ignore:
Timestamp:
2020-03-26T11:57:02+01:00 (2 weeks ago)
Author:
dcarneir
Message:

Inclusion of sea ice thickness in OBS branch

Location:
branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM
Files:
7 edited
1 copied

Legend:

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

    r12140 r12610  
    11981198   ln_sla     = .false.             ! Logical switch for SLA observations 
    11991199   ln_sst     = .false.             ! Logical switch for SST observations 
    1200    ln_sic     = .false.             ! Logical switch for Sea Ice observations 
     1200   ln_sic     = .false.             ! Logical switch for Sea Ice Concentration observations 
     1201   ln_sit     = .false.             ! Logical switch for Sea Ice Thickness observations 
    12011202   ln_vel3d   = .false.             ! Logical switch for velocity observations 
    12021203   ln_sss     = .false.             ! Logical swithc for SSS observations 
     
    12401241   ln_sss_fp_indegs = .true. 
    12411242   ln_sic_fp_indegs = .true. 
     1243   ln_sit_fp_indegs = .true. 
    12421244! All of the *files* variables below are arrays. Use namelist_cfg to add more files 
    12431245   cn_profbfiles = 'profiles_01.nc'      ! Profile feedback input observation file names 
     
    12451247   cn_sstfbfiles = 'sst_01.nc'           ! SST feedback input observation file names 
    12461248   cn_sicfbfiles = 'sic_01.nc'           ! SIC feedback input observation file names 
     1249   cn_sitfbfiles = 'sit_01.nc'           ! SIT feedback input observation file names 
    12471250   cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names 
    12481251   cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names 
     
    12851288   rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees) 
    12861289   rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees) 
     1290   rn_sit_avglamscl = 0.                 ! E/W diameter of SIT observation footprint (metres/degrees) 
     1291   rn_sit_avgphiscl = 0.                 ! N/S diameter of SIT observation footprint (metres/degrees) 
    12871292   nn_1dint = 0                          ! Type of vertical interpolation method 
    12881293   nn_2dint_default = 0                  ! Default horizontal interpolation method 
     
    12911296   nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS 
    12921297   nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC 
     1298   nn_2dint_sit = -1                     ! Horizontal interpolation method for SIT 
    12931299   nn_msshc = 0                          ! MSSH correction scheme 
    12941300   rn_mdtcorr = 1.61                     ! MDT  correction 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r9486 r12610  
    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/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r12140 r12610  
    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 
     
    5354   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
    5455   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
    55    LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     56   LOGICAL :: ln_sic_fp_indegs     !: T=> SIC obs footprint size specified in degrees, F=> in metres 
     57   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
    5658   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology 
    5759   LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 
     
    6567   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
    6668   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
    67    REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
    68    REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
     69   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of SIC observation footprint 
     70   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of SIC observation footprint 
     71   REAL(wp) :: rn_sit_avglamscl     !: E/W diameter of SIT observation footprint 
     72   REAL(wp) :: rn_sit_avgphiscl     !: N/S diameter of SIT observation footprint 
    6973   REAL(wp), PUBLIC :: & 
    7074      &        MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 
     
    7680   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
    7781   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
    78    INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default) 
     82   INTEGER :: nn_2dint_sic     !: SIC horizontal interpolation method (-1 = default) 
     83   INTEGER :: nn_2dint_sit     !: SIT horizontal interpolation method (-1 = default) 
    7984  
    8085   INTEGER, DIMENSION(imaxavtypes) :: & 
     
    136141      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    137142      !!---------------------------------------------------------------------- 
     143#if defined key_cice 
     144USE sbc_oce, ONLY : &        ! CICE variables 
     145   & thick_s                 ! snow depth for freeboard conversion  
     146#endif 
    138147 
    139148      IMPLICIT NONE 
     
    157166         & cn_slafbfiles,      & ! Sea level anomaly input filenames 
    158167         & cn_sicfbfiles,      & ! Seaice concentration input filenames 
     168         & cn_sitfbfiles,      & ! Seaice thickness input filenames 
    159169         & cn_velfbfiles,      & ! Velocity profile input filenames 
    160170         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
     
    189199 
    190200 
     201      LOGICAL :: ln_seaicetypes = .FALSE.          ! Logical switch indicating data type is sea ice 
    191202      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
    192203      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
     
    194205      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
    195206      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     207      LOGICAL :: ln_sit          ! Logical switch for sea ice thickness 
    196208      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
    197209      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     
    252264 
    253265      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    254          &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     266         &            ln_sst, ln_sic, ln_sit, ln_sss, ln_vel3d,               & 
    255267         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
    256268         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
     
    269281         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    270282         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     283         &            ln_sit_fp_indegs,                               & 
    271284         &            cn_profbfiles, cn_slafbfiles,                   & 
    272285         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
     286         &            cn_sitfbfiles,                                  & 
    273287         &            cn_velfbfiles, cn_sssfbfiles,                   & 
    274288         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
     
    292306         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
    293307         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     308         &            rn_sit_avglamscl, rn_sit_avgphiscl,             & 
    294309         &            nn_1dint, nn_2dint_default,                     & 
    295310         &            nn_2dint_sla, nn_2dint_sst,                     & 
    296          &            nn_2dint_sss, nn_2dint_sic,                     & 
     311         &            nn_2dint_sss, nn_2dint_sic, nn_2dint_sit,       & 
    297312         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    298313         &            nn_profdavtypes 
     
    307322      cn_sstfbfiles(:)      = '' 
    308323      cn_sicfbfiles(:)      = '' 
     324      cn_sitfbfiles(:)      = '' 
    309325      cn_velfbfiles(:)      = '' 
    310326      cn_sssfbfiles(:)      = '' 
     
    369385         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
    370386         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
    371          WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
     387         WRITE(numout,*) '             Logical switch for SIC observations                      ln_sic = ', ln_sic 
     388         WRITE(numout,*) '             Logical switch for SIT observations                      ln_sit = ', ln_sit 
    372389         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    373390         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     
    408425         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
    409426         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     427         WRITE(numout,*) '             Type of horizontal interpolation method for SIT    nn_2dint_sit = ', nn_2dint_sit 
    410428         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
    411429         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     
    420438         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
    421439         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
     440         WRITE(numout,*) '             SIT E/W diameter of obs footprint              rn_sit_avglamscl = ', rn_sit_avglamscl 
     441         WRITE(numout,*) '             SIT N/S diameter of obs footprint              rn_sit_avgphiscl = ', rn_sit_avgphiscl 
     442         WRITE(numout,*) '             SIT obs footprint in deg [T] or m [F]          ln_sit_fp_indegs = ', ln_sit_fp_indegs 
    422443         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
    423444         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     
    441462         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
    442463         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
    443       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
     464      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_sss,             & 
    444465         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
    445466         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
     
    560581            cobstypessurf(jtype) = 'sic' 
    561582            clsurffiles(jtype,:) = cn_sicfbfiles 
     583         ENDIF 
     584         IF (ln_sit) THEN 
     585            jtype = jtype + 1 
     586            cobstypessurf(jtype) = 'sit' 
     587            clsurffiles(jtype,:) = cn_sitfbfiles 
    562588         ENDIF 
    563589         IF (ln_sss) THEN 
     
    675701               ztype_avgphiscl = rn_sic_avgphiscl 
    676702               ltype_fp_indegs = ln_sic_fp_indegs 
     703               ltype_night     = .FALSE. 
     704            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
     705               IF ( nn_2dint_sit == -1 ) THEN 
     706                  n2dint_type  = nn_2dint_default 
     707               ELSE 
     708                  n2dint_type  = nn_2dint_sit 
     709               ENDIF 
     710               ztype_avglamscl = rn_sit_avglamscl 
     711               ztype_avgphiscl = rn_sit_avgphiscl 
     712               ltype_fp_indegs = ln_sit_fp_indegs 
    677713               ltype_night     = .FALSE. 
    678714            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     
    863899               & ltype_clim = .TRUE. 
    864900 
    865             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     901            IF ( (TRIM(cobstypessurf(jtype)) == 'sla') .OR. & 
     902               & (TRIM(cobstypessurf(jtype)) == 'sit') ) THEN 
    866903               nvarssurf(jtype) = 1 
    867904               nextrsurf(jtype) = 2 
     
    879916            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
    880917               clvars(1) = 'ICECONC' 
     918            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
     919               clvars(1) = 'FBD' 
     920               ln_seaicetypes = .TRUE. 
    881921            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
    882922               clvars(1) = 'SSS' 
     
    920960               &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 
    921961 
    922             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
     962            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject, ln_seaicetypes ) 
    923963 
    924964            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     
    927967                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    928968            ENDIF 
     969 
     970#if defined key_cice 
     971            IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
     972               CALL obs_rea_snowdepth( surfdataqc(jtype), n2dintsurf(jtype), thick_s(:,:) ) 
     973            ENDIF  
     974#endif 
    929975 
    930976            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
     
    9951041#endif 
    9961042#if defined key_cice 
    997       USE sbc_oce, ONLY : fr_i     ! ice fraction 
     1043      USE sbc_oce, ONLY : &        ! CICE variables 
     1044         & fr_i,          &        ! ice fraction 
     1045         & thick_i                 ! ice thickness 
    9981046#endif 
    9991047#if defined key_top 
     
    13291377                        &           'time-step but some obs are valid then.' ) 
    13301378                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
    1331                         &           ' sea-ice obs will be missed' 
     1379                        &           ' sea-ice concentration obs will be missed' 
    13321380                  ENDIF 
    1333                   surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
    1334                      &                        surfdataqc(jtype)%nsstp(1) 
    1335                   CYCLE 
    13361381               ELSE 
    13371382#if defined key_cice 
     
    13401385                  zsurfvar(:,:) = 1._wp - frld(:,:) 
    13411386#else 
    1342                CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     1387               CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', & 
    13431388                  &           ' but no sea-ice model appears to have been defined' ) 
     1389#endif 
     1390               ENDIF 
     1391 
     1392            CASE('sit') 
     1393               IF ( kstp == 0 ) THEN 
     1394                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
     1395                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
     1396                        &           'time-step but some obs are valid then.' ) 
     1397                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
     1398                        &           ' sea-ice thickness obs will be missed and QC flag set to 4' 
     1399                  ENDIF                                  
     1400               ELSE        
     1401#if defined key_cice 
     1402                  zsurfvar(:,:) = thick_i(:,:) 
     1403#elif defined key_lim2 || defined key_lim3 
     1404                  CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' ) 
     1405#else 
     1406                  CALL ctl_stop( ' Trying to run sea-ice thickness observation operator', & 
     1407                     &           ' but no sea-ice model appears to have been defined' ) 
    13441408#endif 
    13451409               ENDIF 
     
    16411705                  &               lfpindegs(jtype), kmeanstp = imeanstp ) 
    16421706 
     1707 
    16431708            ELSE 
    16441709               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     
    16481713                  &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    16491714                  &               lfpindegs(jtype) ) 
     1715            ENDIF 
     1716 
     1717            ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice 
     1718            ! Thickness") 
     1719            IF ( TRIM(surfdataqc(jtype)%cvars(1)) == 'FBD' ) THEN 
     1720                 surfdata(jtype)%cvars(1) = 'SIT' 
    16501721            ENDIF 
    16511722 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r12140 r12610  
    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 
     
    943947            ENDIF          
    944948         ENDIF 
    945  
     949          
     950         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
     951 
     952#if defined key_cice 
     953         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN 
     954            ! 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)) 
     955            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)  
     956            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2) 
     957            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad 
     958            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN 
     959               surfdataqc%nqc(jobs) = 4 
     960            ENDIF            
     961            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016) 
     962            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ & 
     963                                      (rhow - rhoi) 
     964            ! Flag any negative ice thickness values as bad 
     965            IF (surfdataqc%robs(jobs,1) < 0.0) THEN 
     966               surfdataqc%nqc(jobs) = 4 
     967            ENDIF                                      
     968         ENDIF 
     969          
     970         IF ( zext(1) == obfillflt ) THEN 
     971            ! If the observation value is a fill value, set QC flag to bad 
     972            surfdataqc%nqc(jobs) = 4 
     973         ENDIF 
     974#endif defined key_cice 
     975          
    946976      END DO 
    947977 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r12506 r12610  
    5454 
    5555   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
    56                             kqc_cutoff ) 
     56                            ld_seaicetypes, kqc_cutoff ) 
    5757      !!---------------------------------------------------------------------- 
    5858      !!                    ***  ROUTINE obs_pre_sla  *** 
     
    8484      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
    8585      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
    86       INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
     86      LOGICAL, INTENT(IN) :: ld_seaicetypes        ! Switch to indicate sea ice data 
     87      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff  ! cut off for QC value 
    8788      !! * Local declarations 
    8889      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
     
    141142      ! ----------------------------------------------------------------------- 
    142143 
    143       CALL obs_coo_tim( icycle, & 
    144          &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
    145          &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
    146          &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
    147          &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
     144      IF ( ld_seaicetypes ) THEN 
     145         CALL obs_coo_tim( icycle, & 
     146            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     147            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
     148            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
     149            &              surfdata%nqc,     surfdata%mstp, iotdobs,       & 
     150            &              ld_seaicetypes = ld_seaicetypes ) 
     151      ELSE 
     152         CALL obs_coo_tim( icycle, & 
     153            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     154            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
     155            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
     156            &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
     157      ENDIF 
    148158 
    149159      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     
    559569      &                    kobsno,                                        & 
    560570      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    561       &                    kobsqc,  kobsstp, kotdobs                      ) 
     571      &                    kobsqc,  kobsstp, kotdobs, ld_seaicetypes      ) 
    562572      !!---------------------------------------------------------------------- 
    563573      !!                    ***  ROUTINE obs_coo_tim *** 
     
    607617         & kobsstp          ! Number of time steps up to the  
    608618                            ! observation time 
     619      LOGICAL, OPTIONAL, INTENT(IN) :: ld_seaicetypes 
    609620 
    610621      !! * Local declarations 
     
    621632      INTEGER :: iskip 
    622633      INTEGER :: idaystp 
     634      INTEGER :: icecount 
    623635      REAL(KIND=wp) :: zminstp 
    624636      REAL(KIND=wp) :: zhoustp 
     
    715727            CYCLE 
    716728         ENDIF 
     729 
     730         ! Flag sea ice observations falling on initial timestep 
     731           IF ( PRESENT(ld_seaicetypes) ) THEN 
     732 
     733                IF ( ( kobsstp(jobs) == (nit000 - 1) ) ) THEN 
     734                   IF (lwp) WRITE(numout,*)( 'Sea-ice not initialised on zeroth '// & 
     735                             &    'time-step but SIT observation valid then, flagging '// & 
     736                                  'in time check subroutine obs_coo_tim.' ) 
     737                   kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
     738                   kotdobs      = kotdobs + 1 
     739                   CYCLE 
     740                ENDIF 
     741           ENDIF                      
    717742 
    718743      END DO 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r12140 r12610  
    396396         surfdata%cext(2) = 'MDT' 
    397397      ENDIF 
     398      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'FBD' ) ) THEN 
     399           surfdata%cext(1) = 'freeboard' 
     400           surfdata%cext(2) = 'thick_s' 
     401      ENDIF 
    398402      IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD' 
    399403 
     
    477481               ! Observed value 
    478482               surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
     483               IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN 
     484                   surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
     485                   surfdata%rext(iobs,2) = fbrmdi 
     486               ENDIF 
    479487 
    480488               ! Model and MDT is set to fbrmdi unless read from file 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r11546 r12610  
    530530 
    531531         clfiletype = 'sicfb' 
    532          cllongname = 'Sea ice' 
     532         cllongname = 'Sea ice concentration' 
    533533         clunits    = 'Fraction' 
     534         clgrid     = 'T' 
     535 
     536      CASE('SIT') 
     537 
     538         clfiletype = 'sitfb' 
     539         cllongname = 'Sea ice thickness' 
     540         clunits    = 'm' 
    534541         clgrid     = 'T' 
    535542 
Note: See TracChangeset for help on using the changeset viewer.