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.F90 in NEMO/trunk/src/OCE/ISF – NEMO

source: NEMO/trunk/src/OCE/ISF/isfcav.F90 @ 15058

Last change on this file since 15058 was 15004, checked in by mathiot, 3 years ago

ticket #2960: commit fix to the trunk (WARNING: output convention of isf fluxes changed from oce->isf to isf->oce), no impact on the input file needed for some options
File size: 10.4 KB
Line 
1MODULE isfcav
2   !!======================================================================
3   !!                       ***  MODULE  isfcav  ***
4   !! Ice shelf cavity module :  update ice shelf melting under ice
5   !!                            shelf
6   !!======================================================================
7   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
8   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
9   !!            4.1  !  2019-09  (P. Mathiot) Split ice shelf cavity and ice shelf parametrisation
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   isf_cav       : update ice shelf melting under ice shelf
14   !!----------------------------------------------------------------------
15   USE isf_oce        ! ice shelf public variables
16   !
17   USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
18   USE isfutils , ONLY: debug          ! ice shelf debug subroutine
19   USE isftbl   , ONLY: isf_tbl        ! ice shelf top boundary layer properties subroutine
20   USE isfcavmlt, ONLY: isfcav_mlt     ! ice shelf melt formulation subroutine
21   USE isfcavgam, ONLY: isfcav_gammats ! ice shelf melt exchange coeficient subroutine
22   USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
23   !
24   USE oce      , ONLY: ts                              ! ocean tracers
25   USE par_oce  , ONLY: jpi,jpj                         ! ocean space and time domain
26   USE phycst   , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp  ! physical constants
27   USE eosbn2   , ONLY: ln_teos10                       ! use ln_teos10 or not
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O library
31   USE fldread        ! read input field at current time step
32   USE lbclnk         ! lbclnk
33
34   IMPLICIT NONE
35
36   PRIVATE
37
38   PUBLIC   isf_cav, isf_cav_init ! routine called in isfmlt
39
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE isf_cav( kt, Kmm, ptsc, pqfwf )
48      !!---------------------------------------------------------------------
49      !!                     ***  ROUTINE isf_cav  ***
50      !!
51      !! ** Purpose :   handle surface boundary condition under ice shelf
52      !!
53      !! ** Method  :   based on Mathiot et al. (2017)
54      !!
55      !! ** Action  :   - compute geometry of the Losch top bournary layer (see Losch et al. 2008)
56      !!                - depending on the chooses option
57      !!                   - compute temperature/salt in the tbl
58      !!                   - compute exchange coeficient
59      !!                   - compute heat and fwf fluxes
60      !!                   - output
61      !!
62      !! ** Convention : all fluxes are from isf to oce
63      !!
64      !!---------------------------------------------------------------------
65      !!-------------------------- OUT --------------------------------------
66      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pqfwf  ! ice shelf fwf
67      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc   ! T & S ice shelf cavity contents
68      !!-------------------------- IN  --------------------------------------
69      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
70      INTEGER, INTENT(in) ::   kt   ! ocean time step
71      !!---------------------------------------------------------------------
72      LOGICAL :: lit
73      INTEGER :: nit
74      REAL(wp) :: zerr
75      REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh  ! heat fluxes
76      REAL(wp), DIMENSION(jpi,jpj) :: zqoce_b                  !
77      REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas         ! exchange coeficient
78      REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl             ! temp. and sal. in top boundary layer
79      !!---------------------------------------------------------------------
80      !
81      ! compute T/S/U/V for the top boundary layer
82      CALL isf_tbl(Kmm, ts(:,:,:,jp_tem,Kmm), zttbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
83      CALL isf_tbl(Kmm, ts(:,:,:,jp_sal,Kmm), zstbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
84      !
85      ! output T/S/U/V for the top boundary layer
86      CALL iom_put('ttbl_cav',zttbl(:,:) * mskisf_cav(:,:))
87      CALL iom_put('stbl'    ,zstbl(:,:) * mskisf_cav(:,:))
88      !
89      ! initialisation
90      IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
91      !
92      ! compute ice shelf melting
93      nit = 1 ; lit = .TRUE.
94      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
95         !
96         ! compute gammat everywhere (2d)
97         ! useless if melt specified
98         IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
99            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf,  &
100               &                                    zgammat, zgammas )
101         END IF
102         !
103         ! compute tfrz, latent heat and melt (2d)
104         CALL isfcav_mlt(kt, zgammat, zgammas, zttbl, zstbl, &
105            &                         zqhc   , zqoce, pqfwf  )
106         !
107         ! define if we need to iterate
108         SELECT CASE ( cn_gammablk )
109         CASE ( 'spe','vel' )
110            ! no convergence needed
111            lit = .FALSE.
112         CASE ( 'vel_stab' )
113            ! compute error between 2 iterations
114            zerr = MAXVAL(ABS(zqoce(:,:) - zqoce_b(:,:)))
115            !
116            ! define if iteration needed
117            IF (nit >= 100) THEN              ! too much iteration
118               CALL ctl_stop( 'STOP', 'isf_cav: vel_stab gamma formulation had too many iterations ...' )
119            ELSE IF ( zerr <= 0.01_wp ) THEN  ! convergence is achieve
120               lit = .FALSE.
121            ELSE                              ! converge is not yet achieve
122               nit = nit + 1
123               zqoce_b(:,:) = zqoce(:,:)
124            END IF
125         END SELECT
126         !
127      END DO
128      !
129      ! compute heat and water flux ( > 0 from isf to oce)
130      pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:)
131      zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:)
132      zqhc (:,:) = zqhc(:,:)  * mskisf_cav(:,:)
133      !
134      ! compute heat content flux ( > 0 from isf to oce)
135      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
136      !
137      ! total heat flux ( > 0 from isf to oce)
138      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
139      !
140      ! lbclnk on melt
141      CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp)
142      !
143      ! output fluxes
144      CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
145      !
146      ! set temperature content
147      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
148      !
149      ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
150      IF (lrst_oce) CALL isfrst_write(kt, 'cav', ptsc, pqfwf)
151      !
152      IF ( ln_isfdebug ) THEN
153         IF(lwp) WRITE(numout,*) ''
154         CALL debug('isf_cav: ptsc T',ptsc(:,:,1))
155         CALL debug('isf_cav: ptsc S',ptsc(:,:,2))
156         CALL debug('isf_cav: pqfwf fwf',pqfwf(:,:))
157         IF(lwp) WRITE(numout,*) ''
158      END IF
159      !
160   END SUBROUTINE isf_cav
161
162   SUBROUTINE isf_cav_init
163      !!---------------------------------------------------------------------
164      !!                  ***  ROUTINE isf_cav_init ***
165      !!
166      !! ** Purpose : initialisation of variable needed to compute melt under an ice shelf
167      !!
168      !!----------------------------------------------------------------------
169      INTEGER :: ierr
170      !!---------------------------------------------------------------------
171      !
172      !==============
173      ! 0: allocation
174      !==============
175      !
176      CALL isf_alloc_cav()
177      !
178      !==================
179      ! 1: initialisation
180      !==================
181      !
182      ! top and bottom level of the 'top boundary layer'
183      misfkt_cav(:,:)    = mikt(:,:) ; misfkb_cav(:,:)    = 1
184      !
185      ! thickness of 'tbl' and fraction of bottom cell affected by 'tbl'
186      rhisf_tbl_cav(:,:) = 0.0_wp    ; rfrac_tbl_cav(:,:) = 0.0_wp
187      !
188      ! cavity mask
189      mskisf_cav(:,:)    = (1._wp - tmask(:,:,1)) * ssmask(:,:)
190      !================
191      ! 2: activate restart
192      !================
193      !
194      !================
195      ! 3: read restart
196      !================
197      !
198      ! read cav variable from restart
199      IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
200      !
201      !==========================================
202      ! 3: specific allocation and initialisation (depending of scheme choice)
203      !==========================================
204      !
205      SELECT CASE ( TRIM(cn_isfcav_mlt) )
206      CASE( 'spe' )
207
208         ALLOCATE( sf_isfcav_fwf(1), STAT=ierr )
209         ALLOCATE( sf_isfcav_fwf(1)%fnow(jpi,jpj,1), sf_isfcav_fwf(1)%fdta(jpi,jpj,1,2) )
210         CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_isfdir, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
211
212         IF(lwp) WRITE(numout,*)
213         IF(lwp) WRITE(numout,*) '  ==>> The ice shelf melt inside the cavity is read from forcing files'
214
215      CASE( '2eq' )
216         IF(lwp) WRITE(numout,*)
217         IF(lwp) WRITE(numout,*) '  ==>> The original ISOMIP melt formulation is used to compute melt under the ice shelves'
218
219      CASE( '3eq' )
220         ! coeficient for linearisation of potential tfreez
221         ! Crude approximation for pressure (but commonly used)
222         IF ( ln_teos10 ) THEN   ! linearisation from Jourdain et al. (2017)
223            risf_lamb1 =-0.0564_wp
224            risf_lamb2 = 0.0773_wp
225            risf_lamb3 =-7.8633e-8 * grav * rho0
226         ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
227            risf_lamb1 =-0.0573_wp
228            risf_lamb2 = 0.0832_wp
229            risf_lamb3 =-7.5300e-8 * grav * rho0
230         ENDIF
231
232         IF(lwp) WRITE(numout,*)
233         IF(lwp) WRITE(numout,*) '  ==>> The 3 equations melt formulation is used to compute melt under the ice shelves'
234
235      CASE DEFAULT
236         CALL ctl_stop(' cn_isfcav_mlt method unknown (spe, 2eq, 3eq), check namelist')
237      END SELECT
238      !
239   END SUBROUTINE isf_cav_init
240
241END MODULE isfcav
Note: See TracBrowser for help on using the repository browser.