source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isfcavgam.F90 @ 11889

Last change on this file since 11889 was 11889, checked in by mathiot, 9 months ago

ENHANCE-02_ISF_nemo: add ISOMIP+ test case configuration

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