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/branches/CNRS/dev_r14723_tides_under_isf/src/OCE/ISF – NEMO

source: NEMO/branches/CNRS/dev_r14723_tides_under_isf/src/OCE/ISF/isfcavgam.F90 @ 15605

Last change on this file since 15605 was 15605, checked in by khutchinson, 15 months ago

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