Changeset 13359 for NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90

Timestamp:

2020-07-30T15:40:36+02:00 (4 years ago)

Author:

mathiot

Message:

ticket #1900: add subroutine icb_utl_zavg, fix issue in icb_utl_getkb, add subroutine test_icb_utl_getkb to test awkward cases (will be removed after the review), simplify icbdyn.F90 icbthm.F90 according to the new routine created.

File:

• NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90 (modified) (4 diffs)

NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90

@@ -41 +41 @@
    PUBLIC   icb_utl_copy          ! routine called in icbstp module
    PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules
+   PUBLIC   test_icb_utl_getkb    ! routine called in icbdyn and icbthm modules
+   PUBLIC   icb_utl_zavg          ! routine called in icbdyn and icbthm modules
    PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
    PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module
@@ -482 +484 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(out)               :: kb
+      INTEGER,                INTENT(out):: kb
       REAL(wp), DIMENSION(:), INTENT(in) :: pe3
       REAL(wp),               INTENT(in) :: pD
@@ -488 +490 @@
       INTEGER  :: jk
       REAL(wp) :: zdepw
-      !!
-      zdepw = 0.0
-      kb = 1
+      !!----------------------------------------------------------------------
+      !!
+      zdepw = pe3(1) ; kb = 2
       DO WHILE ( zdepw <  pD)
          zdepw = zdepw + pe3(kb)
          kb = kb + 1
       END DO
-      kb = kb - 1
+      kb = MIN(kb - 1,jpk)
+   END SUBROUTINE
+
+   SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_getkb         ***
+      !!
+      !! ** Purpose :   compute the vertical average of ocean properties affected by icb
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                INTENT(in ) :: kb        ! deepest level affected by icb
+      REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile
+      REAL(wp),               INTENT(in ) :: pD        ! draft
+      REAL(wp),               INTENT(out) :: pzavg     ! z average
+      !!----------------------------------------------------------------------
+      INTEGER  :: jk
+      REAL(wp) :: zdep
+      !!----------------------------------------------------------------------
+      pzavg = 0.0 ; zdep = 0.0
+      DO jk = 1,kb-1
+         pzavg = pzavg + pe3(jk)*pdat(jk)
+         zdep  = zdep  + pe3(jk)
+      END DO
+      ! if kb is limited by mbkt  => bottom value is used between bathy and icb tail
+      ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v)
+      pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD
    END SUBROUTINE
 
@@ -855 +882 @@
    END FUNCTION icb_utl_heat
 
+   SUBROUTINE test_icb_utl_getkb
+      INTEGER :: ikb
+      REAL(wp) :: zD, zout
+      REAL(wp), DIMENSION(jpk) :: ze3, zin
+      WRITE(numout,*) 'Test icb_utl_getkb : '
+      zD = 0.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 8000000.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 80.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 85.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 75.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      WRITE(numout,*) '=================================='
+      WRITE(numout,*) 'Test icb_utl_zavg'
+      zD = 0.0 ; ze3= 20.0 ; zin=1.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+      CALL FLUSH(numout)
+
+      zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      ikb = 2
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      CALL FLUSH(numout)
+
+   END SUBROUTINE test_icb_utl_getkb
+
    !!======================================================================
 END MODULE icbutl