New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
isfcavgam.F90 in NEMO/trunk/src/OCE/ISF – NEMO

source: NEMO/trunk/src/OCE/ISF/isfcavgam.F90 @ 15082

Last change on this file since 15082 was 15082, checked in by clem, 3 years ago

trunk: 1st try to solve problems with ISF option cn_gammablk=vel_stab (partly associated with ticket #2706)
File size: 13.0 KB
RevLine 
[11423]1MODULE isfcavgam
[11395]2   !!======================================================================
3   !!                       ***  MODULE  isfgammats  ***
[11541]4   !! Ice shelf gamma module :  compute exchange coeficient at the ice/ocean interface
[11395]5   !!======================================================================
6   !! History :  4.1  !  (P. Mathiot) original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
[11541]10   !!   isfcav_gammats       : compute exchange coeficient gamma
[11395]11   !!----------------------------------------------------------------------
[12077]12   USE isf_oce
[11852]13   USE isfutils, ONLY: debug
14   USE isftbl  , ONLY: isf_tbl
15
[12068]16   USE oce     , ONLY: uu, vv, rn2         ! ocean dynamics and tracers
[11852]17   USE phycst  , ONLY: grav, vkarmn        ! physical constant
18   USE eosbn2  , ONLY: eos_rab             ! equation of state
19   USE zdfdrg  , ONLY: rCd0_top, r_ke0_top ! vertical physics: top/bottom drag coef.
20   USE iom     , ONLY: iom_put             !
21   USE lib_mpp , ONLY: ctl_stop
22
[11395]23   USE dom_oce        ! ocean space and time domain
24   USE in_out_manager ! I/O manager
25   !
26   IMPLICIT NONE
27   !
28   PRIVATE
29   !
30   PUBLIC   isfcav_gammats
31
[15082]32   !! * Substitutions   
33#  include "do_loop_substitute.h90"
[13237]34#  include "domzgr_substitute.h90"
[11395]35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41   !
42   !!-----------------------------------------------------------------------------------------------------
43   !!                              PUBLIC SUBROUTINES
44   !!-----------------------------------------------------------------------------------------------------
45   !
[12068]46   SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
[11395]47      !!----------------------------------------------------------------------
48      !! ** Purpose    : compute the coefficient echange for heat and fwf flux
49      !!
50      !! ** Method     : select the gamma formulation
[12077]51      !!                 3 method available (cst, vel and vel_stab)
[11395]52      !!---------------------------------------------------------------------
53      !!-------------------------- OUT -------------------------------------
[11852]54      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt  , pgs      ! gamma t and gamma s
[11395]55      !!-------------------------- IN  -------------------------------------
[12068]56      INTEGER                                     :: Kmm             ! ocean time level index
[11395]57      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer
59      !!---------------------------------------------------------------------
60      REAL(wp), DIMENSION(jpi,jpj)                :: zutbl, zvtbl    ! top boundary layer velocity
61      !!---------------------------------------------------------------------
62      !
[11931]63      !==========================================
64      ! 1.: compute velocity in the tbl if needed
65      !==========================================
66      !
[11395]67      SELECT CASE ( cn_gammablk )
[11403]68      CASE ( 'spe'  ) 
69         ! gamma is constant (specified in namelist)
70         ! nothing to do
[12077]71      CASE ('vel', 'vel_stab')
[11495]72         ! compute velocity in tbl
[12068]73         CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
74         CALL isf_tbl(Kmm, vv(:,:,:,Kmm) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav)
[11495]75         !
76         ! mask velocity in tbl with ice shelf mask
[11521]77         zutbl(:,:) = zutbl(:,:) * mskisf_cav(:,:)
78         zvtbl(:,:) = zvtbl(:,:) * mskisf_cav(:,:)
[11495]79         !
80         ! output
[11395]81         CALL iom_put('utbl',zutbl(:,:))
82         CALL iom_put('vtbl',zvtbl(:,:))
83      CASE DEFAULT
84         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
85      END SELECT
86      !
[11931]87      !==========================================
88      ! 2.: compute gamma
89      !==========================================
90      !
[11395]91      SELECT CASE ( cn_gammablk )
92      CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
93         pgt(:,:) = rn_gammat0
94         pgs(:,:) = rn_gammas0
[12077]95      CASE ( 'vel' ) ! gamma is proportional to u*
96         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
97      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
98         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
[11395]99      CASE DEFAULT
100         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
101      END SELECT
102      !
[11931]103      !==========================================
104      ! 3.: output and debug
105      !==========================================
106      !
[11395]107      CALL iom_put('isfgammat', pgt(:,:))
108      CALL iom_put('isfgammas', pgs(:,:))
109      !
[11844]110      IF (ln_isfdebug) THEN
111         CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
112         CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
113      END IF
114      !
[11395]115   END SUBROUTINE isfcav_gammats
116   !
117   !!-----------------------------------------------------------------------------------------------------
118   !!                              PRIVATE SUBROUTINES
119   !!-----------------------------------------------------------------------------------------------------
120   !
[12077]121   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in
[11395]122      &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
123      !!----------------------------------------------------------------------
124      !! ** Purpose    : compute the coefficient echange coefficient
125      !!
[11403]126      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
[11395]127      !!
[12077]128      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
[11395]129      !!---------------------------------------------------------------------
130      !!-------------------------- OUT -------------------------------------
[11852]131      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas [m/s]
[11395]132      !!-------------------------- IN  -------------------------------------
133      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl ! velocity in the losch top boundary layer
134      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd          ! drag coefficient
135      REAL(wp),                     INTENT(in   ) :: pke2         ! background velocity
136      !!---------------------------------------------------------------------
[15082]137      INTEGER  :: ji, jj                     ! loop index
[11395]138      REAL(wp), DIMENSION(jpi,jpj) :: zustar
139      !!---------------------------------------------------------------------
140      !
[15082]141      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
142         ! compute ustar (AD15 eq. 27)
143         zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj)
144         !
145         ! Compute gammats
146         pgt(ji,jj) = zustar(ji,jj) * rn_gammat0
147         pgs(ji,jj) = zustar(ji,jj) * rn_gammas0
148      END_2D
[11395]149      !
150      ! output ustar
151      CALL iom_put('isfustar',zustar(:,:))
152      !
[12077]153   END SUBROUTINE gammats_vel
[11395]154
[12077]155   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
156      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s]
[11395]157      !!----------------------------------------------------------------------
158      !! ** Purpose    : compute the coefficient echange coefficient
159      !!
160      !! ** Method     : gamma is velocity dependent and stability dependent
161      !!
162      !! ** Reference  : Holland and Jenkins, 1999, JPO, p1787-1800
163      !!---------------------------------------------------------------------
164      !!-------------------------- OUT -------------------------------------
[11852]165      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas
[11395]166      !!-------------------------- IN  -------------------------------------
[12068]167      INTEGER                                     :: Kmm            ! ocean time level index
[11395]168      REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared
169      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux
170      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd            ! drag coeficient
171      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl   ! velocity in the losch top boundary layer
172      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl   ! tracer   in the losch top boundary layer
173      !!---------------------------------------------------------------------
174      INTEGER  :: ji, jj                     ! loop index
175      INTEGER  :: ikt                        ! local integer
176      REAL(wp) :: zdku, zdkv                 ! U, V shear
177      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
178      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
179      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
180      REAL(wp) :: zhmax                      ! limitation of mol
181      REAL(wp) :: zetastar                   ! stability parameter
182      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
183      REAL(wp) :: zcoef                      ! temporary coef
184      REAL(wp) :: zdep
185      REAL(wp) :: zeps = 1.0e-20_wp   
186      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
187      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
188      REAL(wp), DIMENSION(2) :: zts, zab
189      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
190      !!---------------------------------------------------------------------
191      !
192      ! compute ustar
[15082]193      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
194         zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )
195      END_2D
[11395]196      !
197      ! output ustar
198      CALL iom_put('isfustar',zustar(:,:))
199      !
200      ! compute Pr and Sc number (eq ??)
201      zPr =   13.8_wp
202      zSc = 2432.0_wp
203      !
204      ! compute gamma mole (eq ??)
205      zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
206      zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
207      !
208      ! compute gamma
[15082]209      DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
210         ikt = mikt(ji,jj)
[11395]211
[15082]212         IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
213            pgt(ji,jj) = rn_gammat0
214            pgs(ji,jj) = rn_gammas0
215         ELSE
216            ! compute Rc number (as done in zdfric.F90)
[11395]217!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
[15082]218            zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
219            !                                            ! shear of horizontal velocity
220            zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
221               &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
222            zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
223               &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
224            !                                            ! richardson number (minimum value set to zero)
225            zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
226           
227            ! compute bouyancy
228            zts(jp_tem) = pttbl(ji,jj)
229            zts(jp_sal) = pstbl(ji,jj)
230            zdep        = gdepw(ji,jj,ikt,Kmm)
231            !
232            CALL eos_rab( zts, zdep, zab, Kmm )
233            !
234            ! compute length scale (Eq ??)
235            zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
236            !
237            ! compute Monin Obukov Length
238            ! Maximum boundary layer depth (Eq ??)
239            zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
240            !
241            ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
242            zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
243            zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
244            !
245            ! compute eta* (stability parameter) (Eq ??)
246            zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
247            !
248            ! compute the sublayer thickness (Eq ??)
249            zhnu = 5 * znu / zustar(ji,jj)
250            !
251            ! compute gamma turb (Eq ??)
252            zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
[11395]253               &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
[15082]254            !
255            ! compute gammats
256            pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
257            pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
258         END IF
259      END_2D
[11395]260
[12077]261   END SUBROUTINE gammats_vel_stab
[11395]262
[11423]263END MODULE isfcavgam
Note: See TracBrowser for help on using the repository browser.