Changeset 11932 for branches

Timestamp:

2019-11-19T18:45:16+01:00 (4 years ago)

Author:

dcarneir

Message:

Changes in OBS and SBC routines for sea ice thickness data assimilation

Location:

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC

Files:

• OBS/obs_oper.F90 (modified) (2 diffs)
• OBS/obs_prep.F90 (modified) (7 diffs)
• OBS/obs_read_surf.F90 (modified) (2 diffs)
• OBS/obs_write.F90 (modified) (1 diff)
• SBC/sbc_ice.F90 (modified) (2 diffs)

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

@@ -37 +37 @@
    USE obs_grid,      ONLY : & 
       & obs_level_search     
+#if defined key_cice
+   USE ice_constants, ONLY : &    ! For conversion from sea ice freeboard to thickness
+      & rhos, rhoi, rhow
+#endif
 
    IMPLICIT NONE
@@ -874 +878 @@
          
          IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
+
+         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
+            ! 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))
+            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1) 
+            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
+            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
+            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
+               surfdataqc%nqc(jobs) = 4
+            ENDIF           
+            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
+            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ &
+                                      (rhow - rhoi)
+            ! Flag any negative ice thickness values as bad
+            IF (surfdataqc%robs(jobs,1) < 0.0) THEN
+               surfdataqc%nqc(jobs) = 4
+            ENDIF                                     
+         ENDIF
          
          IF ( zext(1) == obfillflt ) THEN

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

@@ -53 +53 @@
 
    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
-                            kqc_cutoff )
+                            ld_seaicetypes, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pre_sla  ***
@@ -83 +83 @@
       LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land
       LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary
-      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
+      LOGICAL, INTENT(IN) :: ld_seaicetypes        ! Switch to indicate sea ice data
+      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff  ! cut off for QC value
       !! * Local declarations
       INTEGER :: iqc_cutoff = 255   ! cut off for QC value
@@ -140 +141 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_tim( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &
-         &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &
-         &              surfdata%nqc,     surfdata%mstp, iotdobs        )
+      IF ( ld_seaicetypes ) THEN
+         CALL obs_coo_tim( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &
+            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &
+            &              surfdata%nqc,     surfdata%mstp, iotdobs,       &
+            &              ld_seaicetypes = ld_seaicetypes )
+      ELSE
+         CALL obs_coo_tim( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &
+            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &
+            &              surfdata%nqc,     surfdata%mstp, iotdobs        )
+      ENDIF
 
       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
@@ -558 +568 @@
       &                    kobsno,                                        &
       &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
-      &                    kobsqc,  kobsstp, kotdobs                      )
+      &                    kobsqc,  kobsstp, kotdobs, ld_seaicetypes      )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_tim ***
@@ -606 +616 @@
          & kobsstp          ! Number of time steps up to the 
                             ! observation time
+      LOGICAL, OPTIONAL, INTENT(IN) :: ld_seaicetypes
 
       !! * Local declarations
@@ -620 +631 @@
       INTEGER :: iskip
       INTEGER :: idaystp
+      INTEGER :: icecount
       REAL(KIND=wp) :: zminstp
       REAL(KIND=wp) :: zhoustp
@@ -714 +726 @@
             CYCLE
          ENDIF
+
+         ! Flag sea ice observations falling on initial timestep
+           IF ( PRESENT(ld_seaicetypes) ) THEN
+
+                IF ( ( kobsstp(jobs) == (nit000 - 1) ) ) THEN
+                   IF (lwp) WRITE(numout,*)( 'Sea-ice not initialised on zeroth '// &
+                             &    'time-step but SIT observation valid then, flagging '// &
+                                  'in time check subroutine obs_coo_tim.' )
+                   kobsqc(jobs) = IBSET(kobsqc(jobs),13)
+                   kotdobs      = kotdobs + 1
+                   CYCLE
+                ENDIF
+           ENDIF                     
 
       END DO

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

@@ -380 +380 @@
          surfdata%cext(2) = 'MDT'
       ENDIF
+      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'FBD' ) ) THEN
+           surfdata%cext(1) = 'freeboard'
+           surfdata%cext(2) = 'thick_s'
+      ENDIF
       IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD'
 
@@ -461 +465 @@
                ! Observed value
                surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+               IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN
+                   surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+                   surfdata%rext(iobs,2) = fbrmdi
+               ENDIF
 
                ! Model and MDT is set to fbrmdi unless read from file

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

@@ -530 +530 @@
 
          clfiletype = 'sicfb'
-         cllongname = 'Sea ice'
+         cllongname = 'Sea ice concentration'
          clunits    = 'Fraction'
+         clgrid     = 'T'
+
+      CASE('SIT')
+
+         clfiletype = 'sitfb'
+         cllongname = 'Sea ice thickness'
+         clunits    = 'm'
          clgrid     = 'T'
 

branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -101 +101 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_iu              !: ice thickness at NEMO U point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_iv              !: ice thickness at NEMO V point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_su              !: snow depth at NEMO U point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_sv              !: snow depth at NEMO V point
    
    ! variables used in the coupled interface
@@ -156 +160 @@
                 wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , &
                 ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , &
+                thick_iu(jpi,jpj)     , thick_iv(jpi,jpj)     ,                         &
+                thick_su(jpi,jpj)     , thick_sv(jpi,jpj)     ,                         &
+                ht_i(jpi,jpj,ncat)    , ht_s(jpi,jpj,ncat)    ,                         &
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
                 STAT= ierr(1) )