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

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

Last change on this file since 11844 was 11844, checked in by mathiot, 4 years ago

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