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, 3 years 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.