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/UKMO_MERGE_2019/src/OCE/ISF – NEMO

source: NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavmlt.F90 @ 12199

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

include ENHANCE-02_ISF_nemo in UKMO merge branch

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