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.
isfgammats.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfgammats.F90 @ 11403

Last change on this file since 11403 was 11403, checked in by mathiot, 5 years ago

ENHANCE-02_ISF_nemo : add comments, renaming file (AGRIF), add isfload module (ticket #2142)

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