source: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/isfcavmlt.F90

Last change on this file was 6601, checked in by cetlod, 9 months ago

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