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
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   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
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
32   !! * Substitutions   
33#  include "do_loop_substitute.h90"
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         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 )
99      CASE DEFAULT
100         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
101      END SELECT
102      !
103      !==========================================
104      ! 3.: output and debug
105      !==========================================
106      !
107      CALL iom_put('isfgammat', pgt(:,:))
108      CALL iom_put('isfgammas', pgs(:,:))
109      !
110      IF (ln_isfdebug) THEN
111         CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
112         CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
113      END IF
114      !
115   END SUBROUTINE isfcav_gammats
116   !
117   !!-----------------------------------------------------------------------------------------------------
118   !!                              PRIVATE SUBROUTINES
119   !!-----------------------------------------------------------------------------------------------------
120   !
121   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in
122      &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
123      !!----------------------------------------------------------------------
124      !! ** Purpose    : compute the coefficient echange coefficient
125      !!
126      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
127      !!
128      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
129      !!---------------------------------------------------------------------
130      !!-------------------------- OUT -------------------------------------
131      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas [m/s]
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      !!---------------------------------------------------------------------
137      INTEGER  :: ji, jj                     ! loop index
138      REAL(wp), DIMENSION(jpi,jpj) :: zustar
139      !!---------------------------------------------------------------------
140      !
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
149      !
150      ! output ustar
151      CALL iom_put('isfustar',zustar(:,:))
152      !
153   END SUBROUTINE gammats_vel
154
155   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
156      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s]
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 -------------------------------------
165      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas
166      !!-------------------------- IN  -------------------------------------
167      INTEGER                                     :: Kmm            ! ocean time level index
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
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
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
209      DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
210         ikt = mikt(ji,jj)
211
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)
217!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
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 )) &
253               &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
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
260
261   END SUBROUTINE gammats_vel_stab
262
263END MODULE isfcavgam
Note: See TracBrowser for help on using the repository browser.