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

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

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

ENHANCE-02_ISF_nemo : Initial commit isf simplification (add ISF directory, moved isf routine in and split isf cavity and isf parametrisation, ...) (ticket #2142)

File size: 10.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
15   USE isfdiags
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE eosbn2         ! equation of state
19   USE zdfdrg         ! vertical physics: top/bottom drag coef.
20   !
21   USE in_out_manager ! I/O manager
22   USE iom            ! I/O library
23   USE fldread        ! read input field at current time step
24   USE lbclnk         !
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   isfcav_mlt
30
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
33   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38! -------------------------------------------------------------------------------------------------------
39! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
40! -------------------------------------------------------------------------------------------------------
41
42   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl,   &
43      &                           pqhc, pqoce, pqfwf    )
44      !!----------------------------------------------------------------------
45      !! ** Purpose    :
46      !!
47      !! ** Method     :
48      !!---------------------------------------------------------------------
49      !!-------------------------- OUT -------------------------------------
50      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
51      !!-------------------------- IN  -------------------------------------
52      INTEGER, INTENT(in) :: kt
53      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
54      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
55      !!---------------------------------------------------------------------
56      !!---------------------------------------------------------------------
57      !
58      ! compute latent heat and melt (2d)
59      SELECT CASE ( cn_isfcav_mlt )
60      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
61         CALL isfmlt_spe( kt, pstbl,                  &
62            &                  pqhc, pqoce, pqfwf    )
63      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
64         CALL isfmlt_2eq( pgt, pttbl, pstbl,          &
65            &                  pqhc , pqoce, pqfwf    )
66      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
67         CALL isfmlt_3eq( pgt, pgs , pttbl, pstbl,   &
68            &                  pqhc, pqoce, pqfwf    )
69      CASE ( 'oasis' ) ! fwf pass trough oasis
70         !CALL isfmlt_oasis( kt, pstbl,               &
71         !   &                    zqhc, zqoce, pqfwf  )
72      CASE DEFAULT
73         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
74      END SELECT
75      !
76   END SUBROUTINE isfcav_mlt
77
78! -------------------------------------------------------------------------------------------------------
79! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
80! -------------------------------------------------------------------------------------------------------
81
82   SUBROUTINE isfmlt_spe(kt, pstbl,                  &  ! <<== in
83      &                       pqhc, pqoce, pqfwf    )  ! ==>> out
84      !!----------------------------------------------------------------------
85      !! ** Purpose    :
86      !!
87      !! ** Method     :
88      !!---------------------------------------------------------------------
89      !!-------------------------- OUT -------------------------------------
90      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
91      !!-------------------------- IN  -------------------------------------
92      INTEGER                     , INTENT(in   ) :: kt
93      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl         ! salinity in top boundary layer
94      !!--------------------------------------------------------------------
95      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
96      !!--------------------------------------------------------------------
97      !
98      ! Calculate freezing temperature
99      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
100      !
101      ! read input file
102      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
103      !
104      ! define fwf and qoce
105      pqfwf(:,:) = -sf_isfcav_fwf(1)%fnow(:,:,1)            ! fwf
106      ! ocean heat flux is assume to be equal to the latent heat
107      pqoce(:,:) =  pqfwf(:,:) * rLfusisf             ! ocean heat flux
108      pqhc (:,:) =  pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux
109      !
110   END SUBROUTINE isfmlt_spe
111
112   SUBROUTINE isfmlt_2eq(pgt, pttbl, pstbl,          &  ! <<== in
113      &                       pqhc , pqoce, pqfwf    )  ! ==>> out
114      !!----------------------------------------------------------------------
115      !! ** Purpose    :
116      !!
117      !! ** Method     :
118      !!---------------------------------------------------------------------
119      !!-------------------------- OUT -------------------------------------
120      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
121      !!-------------------------- IN  -------------------------------------
122      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
123      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
124      !!--------------------------------------------------------------------
125      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
126      !!--------------------------------------------------------------------
127      !
128      ! Calculate freezing temperature
129      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
130      !
131      ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
132      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * ( pttbl(:,:)-ztfrz(:,:) ) * r1_Lfusisf  ! fresh water flux ( > 0 out )
133      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                                         ! ocea-ice flux (assume to be equal to latent heat flux) ( > 0 out )
134      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                                 ! heat content flux ( > 0 out )
135      !
136      ! output thermal driving
137      CALL iom_put('isfthermald_cav',( pttbl(:,:)-ztfrz(:,:) ))
138      !
139   END SUBROUTINE isfmlt_2eq
140
141   SUBROUTINE isfmlt_3eq(pgt, pgs , pttbl, pstbl,   &
142      &                       pqhc, pqoce, pqfwf    )
143      !!----------------------------------------------------------------------
144      !! ** Purpose    :
145      !!
146      !! ** Method     :
147      !!---------------------------------------------------------------------
148      !!-------------------------- OUT -------------------------------------
149      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
150      !!-------------------------- IN  -------------------------------------
151      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! temperature exchange coeficient
152      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
153      !!--------------------------------------------------------------------
154      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
155      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
156      REAL(wp) :: zeps = 1.e-20
157      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz             ! freezing point
158      REAL(wp), DIMENSION(jpi,jpj) :: zqcon             ! conductive flux through the ice shelf
159      !
160      INTEGER  ::   ji, jj     ! dummy loop indices
161      !!--------------------------------------------------------------------
162      !
163      ! compute upward heat flux zhtflx and upward water flux zwflx
164      ! Resolution of a 2d equation from equation 24, 25 and 26
165      DO jj = 1, jpj
166         DO ji = 1, jpi
167            !
168            ! compute coeficient to solve the 2nd order equation
169            zeps1 = rau0_rcp*pgt(ji,jj)
170            zeps2 = rLfusisf*rau0*pgs(ji,jj)
171            zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
172            zeps4 = risf_lamb2+risf_lamb3*risfdep(ji,jj)
173            zeps6 = zeps4-pttbl(ji,jj)
174            zeps7 = zeps4-rtsurf
175            !
176            ! solve the 2nd order equation to find zsfrz
177            zaqe  = risf_lamb1 * (zeps1 + zeps3)
178            zaqer = 0.5_wp/MIN(zaqe,-risf_eps)
179            zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
180            zcqe  = zeps2*pstbl(ji,jj)
181            zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
182            !
183            ! Presumably zdis can never be negative because gammas is very small compared to gammat
184            zsfrz=(-zbqe-SQRT(zdis))*zaqer
185            IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer  ! check this if this if is needed
186            !
187            ! compute t freeze (eq. 25)
188            ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
189            !
190            ! compute the upward water and heat flux (eq. 24 and eq. 26)
191            ! ocean heat content flux added later on.
192            pqfwf(ji,jj) = rau0 * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux (> 0 out)
193            pqoce(ji,jj) = pgt(ji,jj) * rau0_rcp * ( pttbl(ji,jj) - ztfrz(ji,jj) )        ! ocean-ice heat flux (> 0 out)
194            pqhc (ji,jj) = pqfwf(ji,jj) * ztfrz(ji,jj) * rcp                              ! heat content   flux (> 0 out)
195            zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf ) ! to be check
196            !
197         END DO
198      END DO
199      !
200      ! output conductive heat flux through the ice
201      CALL iom_put('qconisf', zqcon)
202      !
203      ! output thermal driving
204      CALL iom_put('isfthermald_cav',( pttbl(:,:) - ztfrz(:,:) ))
205      !
206   END SUBROUTINE isfmlt_3eq
207
208   !SUBROUTINE isfmlt_3eq_frz_ktm1
209   ! compute tfrz based on sfrz value at kt-1 (need to be SAVED local array)
210   ! compute qfwf (eq 24)
211   ! compute zqoce, zqlat, zqcon, zqhc
212   ! compute sfrz (eq 26)
213   !END SUBROUTINE isfmlt_3eq_frz_ktm1
214
215END MODULE isfcavmlt
Note: See TracBrowser for help on using the repository browser.