Changeset 9513

Timestamp:

2018-04-27T14:30:02+02:00 (6 years ago)

Author:

mathiot

Message:

Add option to detect and remove subglacial lake (do not affect closed sea option)

Location:

branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM

Files:

• CONFIG/SHARED/namelist_ref (modified) (1 diff)
• NEMO/OPA_SRC/DOM/dommsk.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/DOM/domngb.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/DOM/domutl.F90 (added)
• NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (2 diffs)

branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/CONFIG/SHARED/namelist_ref

@@ -129 +129 @@
    rn_zisfmax       = 7000.       ! limiter in m     on the previous condition. (if change larger than this number, get back to value before we enforce the connection)
    ln_isfcheminey   = .false.     ! close cheminey
+   ln_isfsubgl      = .true.      ! remove subglacial lake created by the remapping process
+   rn_isfsubgllon   =    0.0      !  longitude of the seed to determine the open ocean
+   rn_isfsubgllat   =    0.0      !  latitude  of the seed to determine the open ocean
 /
 !-----------------------------------------------------------------------

branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

@@ -25 +25 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE domngb          ! find nearest wet point
+   USE domutl          ! fill closed 3d pool below isf
+   USE domzgr, ONLY : ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat ! import ln_isfsubgl to mask close sea below isf
+   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -133 +137 @@
       INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
       INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
+      INTEGER  ::   jiseed, jjseed           !   -       -
       INTEGER  ::   ios
       INTEGER  ::   isrow                    ! index for ORCA1 starting row
@@ -187 +192 @@
       END DO  
       
-      ! (ISF) define barotropic mask and mask the ice shelf point
-      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
-      
       DO jk = 1, jpk
          DO jj = 1, jpj
@@ -199 +201 @@
          END DO  
       END DO  
-
+      !
+      IF ( ln_isfsubgl ) THEN
+         ! check closed wet pool
+         CALL dom_ngb(rn_isfsubgllon, rn_isfsubgllat, jiseed, jjseed, 'T', lwet=.TRUE.)
+         CALL fill_pool( jiseed, jjseed, 1, tmask, -1._wp )
+         ! at this point itab3d (:,1:ijmax,:) can have 3 different values :
+         !              0 where there where already 0
+         !              -1 where the ocean points are connected
+         !              1 where ocean points in tmask are not connected
+         IF (lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*)'dommsk : removal of subglacial lakes '
+            WRITE(numout,*)'~~~~~~~'
+            WRITE(numout,*)'Number of disconected points : ', COUNT(  (tmask(:,:,:) == 1) )
+            WRITE(numout,*)'lon/lat seed to detect main ocean is: ', rn_isfsubgllon, rn_isfsubgllat
+            WRITE(numout,*)'i/j     seed to detect main ocean is: ', jiseed, jjseed
+         END IF
+         DO jk = 1, jpk
+            WHERE (tmask(:,:,jk) > 0 .AND. misfdep(:,:) > 1) tmask(:,:,jk) = 0 ! remove only subglacial lake (ie similar to close sea only below an ice shelf 
+            WHERE (tmask(:,:,jk) < 0) tmask(:,:,jk) = 1 ! restore mask value
+         END DO
+      END IF
+
+      ! (ISF) define barotropic mask and mask the ice shelf point
+      DO jj = 1, jpj
+         DO ji = 1, jpi   ! vector loop
+            ssmask(ji,jj)  = MIN(1._wp,SUM(tmask(ji,jj,:)))
+         END DO
+      END DO
 !!gm  ????
 #if defined key_zdfkpp

branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

@@ -28 +28 @@
 CONTAINS
 
-   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )
+   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, lwet, kkk )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dom_ngb  ***
       !!
-      !! ** Purpose :   find the closest grid point from a given lon/lat position
+      !! ** Purpose :   find the closest wet grid point from a given lon/lat position
       !!
       !! ** Method  :   look for minimum distance in cylindrical projection 
@@ -40 +40 @@
       REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
       INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
+      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
+      LOGICAL         , INTENT(in   ), OPTIONAL :: lwet ! logical to decide if we look for every where (false) 
+                                                                !    or only over the wet point (true)
       CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W'
       !
+      INTEGER :: ik         ! working level
       INTEGER , DIMENSION(2) ::   iloc
       REAL(wp)               ::   zlon, zmini
       REAL(wp), POINTER, DIMENSION(:,:) ::  zglam, zgphi, zmask, zdist
+      LOGICAL :: llwet      ! working logical
       !!--------------------------------------------------------------------
       !
@@ -51 +56 @@
       CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist )
       !
+      ! select lat/lon variable
+      SELECT CASE( cdgrid )
+         CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:)
+         CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:)
+         CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:)
+         CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:)
+      END SELECT
+      !
+      ! select vertical level
+      ik   = 1
+      IF ( PRESENT(kkk) ) ik=kkk
+      !
+      ! select mask variable
       zmask(:,:) = 0._wp
-      SELECT CASE( cdgrid )
-      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1)
-      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1)
-      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1)
-      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1)
-      END SELECT
-
-      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
-      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
-      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
-      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
-
-      zglam(:,:) = zglam(:,:) - zlon
+      llwet = .TRUE.
+      IF ( PRESENT(lwet)) llwet=lwet 
+      IF (llwet) THEN
+         SELECT CASE( cdgrid )
+            CASE( 'U' )  ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
+            CASE( 'V' )  ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
+            CASE( 'F' )  ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
+            CASE DEFAULT ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
+         END SELECT
+      ELSE
+         zmask(nldi:nlei,nldj:nlej)=1.e0
+      END IF
+      !
+      ! compute distance
+      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN
+         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
+         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
+         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
+         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
+         zglam(:,:) = zglam(:,:) - zlon
+      ELSE
+         zglam(:,:) = zglam(:,:) - plon
+      END IF
+      !
       zgphi(:,:) = zgphi(:,:) - plat
       zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
-      
+      !
+      ! find min
       IF( lk_mpp ) THEN  
          CALL mpp_minloc( zdist(:,:), zmask, zmini, kii, kjj)

branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -50 +50 @@
    PUBLIC   dom_zgr        ! called by dom_init.F90
 
+   ! remove subglacial lake under the ice sheet.
+   LOGICAL , PUBLIC  :: ln_isfsubgl
+   REAL(wp), PUBLIC  :: rn_isfsubgllon, rn_isfsubgllat
    !                              !!* Namelist namzgr_sco *
    LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
@@ -1293 +1296 @@
       ! NAMELIST
       INTEGER  ::   ios
-      INTEGER  :: nn_kisfmax    = 999. !  
-      REAL(wp) :: rn_isfdep_min = 10.0_wp ! ice shelf minimal thickness 
+      INTEGER  :: nn_kisfmax    = 999.                              !  
+      REAL(wp) :: rn_isfdep_min = 10.0_wp                           ! ice shelf minimal thickness 
       REAL(wp) :: rn_isfhw_deep = 1.e-3 , rn_isfhw_shallow = 1.0e-3 ! threshold to define grounding line in deep/shallow water
-      REAL(wp) :: rn_isfshallow = 0.0_wp  ! threshold to define shallow ice shelf cavity
-      REAL(wp) :: rn_zisfmax    = 6000.0_wp ! maximun meter of ice we are allowed to dig to assure connectivity
-      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity
+      REAL(wp) :: rn_isfshallow = 0.0_wp                            ! threshold to define shallow ice shelf cavity
+      REAL(wp) :: rn_zisfmax    = 6000.0_wp                         ! maximun meter of ice we are allowed to dig to assure connectivity
+      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity, remove channel in water column (on z space not level space)
       !!---------------------------------------------------------------------
       NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
               &           rn_isfdep_min, rn_isfhw_deep, rn_isfhw_shallow, rn_isfshallow, &
-              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel
+              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
       !
       IF( nn_timing == 1 )  CALL timing_start('zgr_isf')