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

Last change on this file since 11423 was 11423, checked in by mathiot, 15 months ago

ENHANCE-02_ISF_nemo : add UKESM ice sheet coupling method (ticket #2142)

File size: 15.9 KB
Line 
1MODULE isfcavmlt
2   !!======================================================================
3   !!                       ***  MODULE  isfcav_mlt  ***
4   !! ice shelf module :  update surface ocean boundary condition under ice
5   !!                   shelf
6   !!======================================================================
7   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   isfcav_mlt    :
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE isf            ! ice shelf public variables
15   USE dom_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   USE eosbn2         ! equation of state
18   !
19   USE in_out_manager ! I/O manager
20   USE iom            ! I/O library
21   USE fldread        ! read input field at current time step
22   USE lib_fortran
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   isfcav_mlt
28
29   !!----------------------------------------------------------------------
30   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
31   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
32   !! Software governed by the CeCILL license (see ./LICENSE)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36! -------------------------------------------------------------------------------------------------------
37! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
38! -------------------------------------------------------------------------------------------------------
39
40   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
41      &                           pqhc, pqoce, pqfwf  )
42      !!----------------------------------------------------------------------
43      !!
44      !! ** Purpose    : compute ice shelf fwf/heqt fluxes
45      !!
46      !!---------------------------------------------------------------------
47      !!-------------------------- OUT -------------------------------------
48      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
49      !!-------------------------- IN  -------------------------------------
50      INTEGER, INTENT(in) :: kt
51      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
52      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
53      !!---------------------------------------------------------------------
54      !!---------------------------------------------------------------------
55      !
56      ! compute latent heat and melt (2d)
57      SELECT CASE ( cn_isfcav_mlt )
58      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
59         CALL isfmlt_spe( kt, pstbl,               &
60            &                  pqhc, pqoce, pqfwf  )
61      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
62         CALL isfmlt_2eq( pgt, pttbl, pstbl,       &
63            &                  pqhc , pqoce, pqfwf )
64      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
65         CALL isfmlt_3eq( pgt, pgs , pttbl, pstbl, &
66            &                  pqhc, pqoce, pqfwf  )
67      CASE ( 'oasis' ) ! fwf pass trough oasis
68         CALL isfmlt_oasis( kt, pstbl,              &
69            &                   pqhc, pqoce, pqfwf  )
70      CASE DEFAULT
71         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
72      END SELECT
73      !
74   END SUBROUTINE isfcav_mlt
75
76! -------------------------------------------------------------------------------------------------------
77! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
78! -------------------------------------------------------------------------------------------------------
79
80   SUBROUTINE isfmlt_spe(kt, pstbl,              &  ! <<== in
81      &                      pqhc , pqoce, pqfwf )  ! ==>> out
82      !!----------------------------------------------------------------------
83      !! ** Purpose    : - read ice shelf melt from forcing file
84      !!                 - compute ocea-ice heat flux (assuming it is equal to latent heat)
85      !!                 - compute heat content flux
86      !!---------------------------------------------------------------------
87      !!-------------------------- OUT -------------------------------------
88      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
89      !!-------------------------- IN  -------------------------------------
90      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
91      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
92      !!--------------------------------------------------------------------
93      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
94      !!--------------------------------------------------------------------
95      !
96      ! Calculate freezing temperature
97      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
98      !
99      ! read input file
100      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
101      !
102      ! define fwf and qoce
103      ! ocean heat flux is assume to be equal to the latent heat
104      pqfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)      ! fwf                ( >0 out)
105      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
106      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
107      !
108   END SUBROUTINE isfmlt_spe
109
110   SUBROUTINE isfmlt_2eq(pgt, pttbl, pstbl,       &  ! <<== in
111      &                       pqhc , pqoce, pqfwf )  ! ==>> out
112      !!----------------------------------------------------------------------
113      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
114      !!
115      !! ** Method     : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
116      !!                 From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
117      !!                   qfwf  = Gammat * Rau0 * Cp * ( Tw - Tfrz ) / Lf
118      !!                   qhoce = qlat
119      !!                   qhc   = qfwf * Cp * Tfrz
120      !!
121      !! ** Reference  : Hunter,  J.  R.:  Specification  for  test  models  of  ice  shelf  cavities, 
122      !!                 Tech.  Rep.  June,  Antarctic  Climate  &  Ecosystems  Cooperative  Research  Centre,  available  at: 
123      !!                 http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
124      !!---------------------------------------------------------------------
125      !!-------------------------- OUT -------------------------------------
126      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
127      !!-------------------------- IN  -------------------------------------
128      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
129      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
130      !!--------------------------------------------------------------------
131      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
132      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
133      !!--------------------------------------------------------------------
134      !
135      ! Calculate freezing temperature
136      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
137      !
138      ! thermal driving
139      zthd (:,:) = pttbl(:,:) - ztfrz(:,:)
140      !
141      ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
142      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) * r1_Lfusisf  ! fresh water flux  ( > 0 out )
143      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
144      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
145      !
146      ! output thermal driving
147      CALL iom_put('isfthermald_cav', zthd )
148      !
149   END SUBROUTINE isfmlt_2eq
150
151   SUBROUTINE isfmlt_3eq(pgt, pgs , pttbl, pstbl, &
152      &                       pqhc, pqoce, pqfwf  )
153      !!----------------------------------------------------------------------
154      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation
155      !!
156      !! ** Method     : The melt rate is determined considering the heat balance, the salt balance
157      !!                 at the phase change interface and a linearisation of the equation of state.
158      !!
159      !! ** Reference  : - Holland, D. M. and Jenkins, A.,
160      !!                   Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf,
161      !!                   J. Phys. Oceanogr., 29, 1999.
162      !!                 - Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone,
163      !!                   R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland,
164      !!                   P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.:
165      !!                   Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects:
166      !!                   MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1),
167      !!                   Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016.
168      !!---------------------------------------------------------------------
169      !!-------------------------- OUT -------------------------------------
170      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
171      !!-------------------------- IN  -------------------------------------
172      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs          ! heat/salt exchange coeficient
173      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl        ! mean temperature and salinity in top boundary layer
174      !!--------------------------------------------------------------------
175      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7       ! dummy local scalar for quadratic equation resolution
176      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac     ! dummy local scalar for quadratic equation resolution
177      REAL(wp) :: zeps = 1.e-20
178      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing point
179      REAL(wp), DIMENSION(jpi,jpj) :: zqcon         ! conductive flux through the ice shelf
180      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
181      !
182      INTEGER  ::   ji, jj     ! dummy loop indices
183      !!--------------------------------------------------------------------
184      !
185      ! compute upward heat flux zhtflx and upward water flux zwflx
186      ! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
187      DO jj = 1, jpj
188         DO ji = 1, jpi
189            !
190            ! compute coeficient to solve the 2nd order equation
191            zeps1 = rau0_rcp*pgt(ji,jj)
192            zeps2 = rLfusisf*rau0*pgs(ji,jj)
193            zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
194            zeps4 = risf_lamb2+risf_lamb3*risfdep(ji,jj)
195            zeps6 = zeps4-pttbl(ji,jj)
196            zeps7 = zeps4-rtsurf
197            !
198            ! solve the 2nd order equation to find zsfrz
199            zaqe  = risf_lamb1 * (zeps1 + zeps3)
200            zaqer = 0.5_wp/MIN(zaqe,-zeps)
201            zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
202            zcqe  = zeps2*pstbl(ji,jj)
203            zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
204            !
205            ! Presumably zdis can never be negative because gammas is very small compared to gammat
206            zsfrz=(-zbqe-SQRT(zdis))*zaqer
207            IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer  ! check this if this if is needed
208            !
209            ! compute t freeze (eq. 25)
210            ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
211            !
212            ! thermal driving
213            zthd(ji,jj) = pttbl(ji,jj) - ztfrz(ji,jj)
214            !
215            ! compute the upward water and heat flux (eq. 24 and eq. 26)
216            pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
217            pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
218            pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out)
219            !
220            zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )                                  ! conductive flux through the ice (> 0 out)
221            !
222         END DO
223      END DO
224      !
225      ! output conductive heat flux through the ice
226      CALL iom_put('qconisf', zqcon)
227      !
228      ! output thermal driving
229      CALL iom_put('isfthermald_cav', zthd)
230      !
231   END SUBROUTINE isfmlt_3eq
232
233   SUBROUTINE isfmlt_oasis(kt, pstbl,              &  ! <<== in
234      &                        pqhc , pqoce, pqfwf )  ! ==>> out
235      !!----------------------------------------------------------------------
236      !!
237      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
238      !!
239      !! ** Purpose    : - read ice shelf melt from forcing file => pattern
240      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
241      !!                 - scale fwf and compute heat fluxes
242      !!
243      !!---------------------------------------------------------------------
244      !!-------------------------- OUT -------------------------------------
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
246      !!-------------------------- IN  -------------------------------------
247      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
248      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
249      !!--------------------------------------------------------------------
250      REAL(wp)                     :: zfwf_fld, zfwf_oasis               ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
251      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
252      REAL(wp), DIMENSION(jpi,jpj) :: zfwf                               ! 2d fwf map after scaling
253      !!--------------------------------------------------------------------
254      !
255      ! Calculate freezing temperature
256      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
257      !
258      ! read input file
259      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
260      !
261      ! ice shelf 2d map
262      zfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)
263      !
264      ! compute glob sum from input file
265      zfwf_fld = glob_sum('isfcav_mlt', zfwf(:,:))
266      !
267      ! compute glob sum from atm->oce ice shelf fwf
268      zfwf_oasis = glob_sum('isfcav_mlt', fwfisf_oasis(:,:))
269      !
270      ! scale fwf
271      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
272      !
273      ! define fwf and qoce
274      ! ocean heat flux is assume to be equal to the latent heat
275      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out)
276      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
277      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
278      !
279   END SUBROUTINE isfmlt_oasis
280
281   !SUBROUTINE isfmlt_3eq_frz_ktm1
282   ! compute tfrz based on sfrz value at kt-1 (need to be SAVED local array)
283   ! => should reduce error due to linarisation
284   ! compute qfwf (eq 24)
285   ! compute zqoce, zqlat, zqcon, zqhc
286   ! compute sfrz (eq 26)
287   !END SUBROUTINE isfmlt_3eq_frz_ktm1
288
289END MODULE isfcavmlt
Note: See TracBrowser for help on using the repository browser.