Changeset 15084 for NEMO/trunk/src/OCE/ISF/isfcav.F90

Timestamp:

2021-07-05T21:48:42+02:00 (3 years ago)

Author:

clem

Message:

trunk ISF: correct option cn_gammablk=vel_stab as much as I understand it and remove some useless lbc_lnk. Ref ticket is #2706

File:

• NEMO/trunk/src/OCE/ISF/isfcav.F90 (modified) (8 diffs)

NEMO/trunk/src/OCE/ISF/isfcav.F90

@@ -22 +22 @@
    USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
    !
-   USE oce      , ONLY: ts                              ! ocean tracers
+   USE oce      , ONLY: ts, uu, vv, rn2                 ! ocean dynamics and tracers
    USE par_oce  , ONLY: jpi,jpj                         ! ocean space and time domain
    USE phycst   , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp  ! physical constants
@@ -31 +31 @@
    USE fldread        ! read input field at current time step
    USE lbclnk         ! lbclnk
+   USE lib_mpp        ! MPP library
 
    IMPLICIT NONE
@@ -38 +39 @@
    PUBLIC   isf_cav, isf_cav_init ! routine called in isfmlt
 
+   !! * Substitutions   
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -71 +74 @@
       !!---------------------------------------------------------------------
       LOGICAL :: lit
-      INTEGER :: nit
+      INTEGER :: nit, ji, jj, ikt
       REAL(wp) :: zerr
+      REAL(wp) :: zcoef, zdku, zdkv
       REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh  ! heat fluxes
-      REAL(wp), DIMENSION(jpi,jpj) :: zqh_b                    !
+      REAL(wp), DIMENSION(jpi,jpj) :: zqh_b, zRc               !
       REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas         ! exchange coeficient
       REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl             ! temp. and sal. in top boundary layer
@@ -91 +95 @@
          zqoce(:,:) = -pqfwf(:,:) * rLfusisf ! 
          zqh_b(:,:) =  ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
+
+         DO_2D( 0, 0, 0, 0 )
+            ikt = mikt(ji,jj)
+            ! compute Rc number (as done in zdfric.F90)
+!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
+            zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
+            !                                            ! shear of horizontal velocity
+            zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
+               &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
+            zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
+               &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
+            !                                            ! richardson number (minimum value set to zero)
+            zRc(ji,jj) = MAX(rn2(ji,jj,ikt+1), 1.e-20_wp) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
+         END_2D
+         CALL lbc_lnk( 'isfmlt', zRc, 'T', 1._wp )
       ENDIF
       !
@@ -100 +119 @@
          ! useless if melt specified
          IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
-            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf,  &
+            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf, zRc,  &
                &                                    zgammat, zgammas )
          END IF
@@ -115 +134 @@
          CASE ( 'vel_stab' )
             ! compute error between 2 iterations
-            zerr = MAXVAL(ABS(zqhc(:,:)+zqoce(:,:) - zqh_b(:,:)))
+            zerr = 0._wp
+            DO_2D( 0, 0, 0, 0 )
+               zerr = MAX( zerr, ABS(zqhc(ji,jj)+zqoce(ji,jj) - zqh_b(ji,jj)) )
+            END_2D
+            CALL mpp_max( 'isfcav', zerr )   ! max over the global domain
             !
             ! define if iteration needed
@@ -130 +153 @@
       END DO
       !
-      ! compute heat and water flux ( > 0 from isf to oce)
-      pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:)
-      zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:)
-      zqhc (:,:) = zqhc(:,:)  * mskisf_cav(:,:)
-      !
-      ! compute heat content flux ( > 0 from isf to oce)
-      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
-      !
-      ! total heat flux ( > 0 from isf to oce)
-      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
-      !
-      ! lbclnk on melt
-      CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ! compute heat and water flux ( > 0 from isf to oce)
+         pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_cav(ji,jj)
+         zqoce(ji,jj) = zqoce(ji,jj) * mskisf_cav(ji,jj)
+         zqhc (ji,jj) = zqhc(ji,jj)  * mskisf_cav(ji,jj)
+         !
+         ! compute heat content flux ( > 0 from isf to oce)
+         zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf    ! 2d latent heat flux (W/m2)
+         !
+         ! total heat flux ( > 0 from isf to oce)
+         zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
+         !
+         ! set temperature content
+         ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
+      END_2D
       !
       ! output fluxes
       CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
-      !
-      ! set temperature content
-      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
       !
       ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)