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 @ 15085

Last change on this file since 15085 was 15085, checked in by clem, 3 years ago

ISF further cleaning and debugging(?) for isfcpl but I think it requires more attention

File size: 11.7 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, uu, vv, rn2                 ! ocean dynamics and 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   USE lib_mpp        ! MPP library
34
35   IMPLICIT NONE
36
37   PRIVATE
38
39   PUBLIC   isf_cav, isf_cav_init ! routine called in isfmlt
40
41   !! * Substitutions   
42#  include "do_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE isf_cav( kt, Kmm, ptsc, pqfwf )
51      !!---------------------------------------------------------------------
52      !!                     ***  ROUTINE isf_cav  ***
53      !!
54      !! ** Purpose :   handle surface boundary condition under ice shelf
55      !!
56      !! ** Method  :   based on Mathiot et al. (2017)
57      !!
58      !! ** Action  :   - compute geometry of the Losch top bournary layer (see Losch et al. 2008)
59      !!                - depending on the chooses option
60      !!                   - compute temperature/salt in the tbl
61      !!                   - compute exchange coeficient
62      !!                   - compute heat and fwf fluxes
63      !!                   - output
64      !!
65      !! ** Convention : all fluxes are from isf to oce
66      !!
67      !!---------------------------------------------------------------------
68      !!-------------------------- OUT --------------------------------------
69      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pqfwf  ! ice shelf fwf
70      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc   ! T & S ice shelf cavity contents
71      !!-------------------------- IN  --------------------------------------
72      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
73      INTEGER, INTENT(in) ::   kt   ! ocean time step
74      !!---------------------------------------------------------------------
75      LOGICAL :: lit
76      INTEGER :: nit, ji, jj, ikt
77      REAL(wp) :: zerr
78      REAL(wp) :: zcoef, zdku, zdkv
79      REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh  ! heat fluxes
80      REAL(wp), DIMENSION(jpi,jpj) :: zqh_b, zRc               !
81      REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas         ! exchange coeficient
82      REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl             ! temp. and sal. in top boundary layer
83      !!---------------------------------------------------------------------
84      !
85      ! compute T/S/U/V for the top boundary layer
86      CALL isf_tbl(Kmm, ts(:,:,:,jp_tem,Kmm), zttbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
87      CALL isf_tbl(Kmm, ts(:,:,:,jp_sal,Kmm), zstbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
88      !
89      ! output T/S/U/V for the top boundary layer
90      CALL iom_put('ttbl_cav',zttbl(:,:) * mskisf_cav(:,:))
91      CALL iom_put('stbl'    ,zstbl(:,:) * mskisf_cav(:,:))
92      !
93      ! initialisation
94      IF ( TRIM(cn_gammablk) == 'vel_stab' ) THEN
95         zqoce(:,:) = -pqfwf(:,:) * rLfusisf !
96         zqh_b(:,:) =  ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
97
98         DO_2D( 0, 0, 0, 0 )
99            ikt = mikt(ji,jj)
100            ! compute Rc number (as done in zdfric.F90)
101!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
102            zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
103            !                                            ! shear of horizontal velocity
104            zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
105               &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
106            zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
107               &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
108            !                                            ! richardson number (minimum value set to zero)
109            zRc(ji,jj) = MAX(rn2(ji,jj,ikt+1), 1.e-20_wp) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
110         END_2D
111         CALL lbc_lnk( 'isfcav', zRc, 'T', 1._wp )
112      ENDIF
113      !
114      ! compute ice shelf melting
115      nit = 1 ; lit = .TRUE.
116      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
117         !
118         ! compute gammat everywhere (2d)
119         ! useless if melt specified
120         IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
121            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf, zRc,  &
122               &                                    zgammat, zgammas )
123         END IF
124         !
125         ! compute tfrz, latent heat and melt (2d)
126         CALL isfcav_mlt(kt, zgammat, zgammas, zttbl, zstbl, &
127            &                         zqhc   , zqoce, pqfwf  )
128         !
129         ! define if we need to iterate
130         SELECT CASE ( cn_gammablk )
131         CASE ( 'spe','vel' )
132            ! no convergence needed
133            lit = .FALSE.
134         CASE ( 'vel_stab' )
135            ! compute error between 2 iterations
136            zerr = 0._wp
137            DO_2D( 0, 0, 0, 0 )
138               zerr = MAX( zerr, ABS(zqhc(ji,jj)+zqoce(ji,jj) - zqh_b(ji,jj)) )
139            END_2D
140            CALL mpp_max( 'isfcav', zerr )   ! max over the global domain
141            !
142            ! define if iteration needed
143            IF (nit >= 100) THEN              ! too much iteration
144               CALL ctl_stop( 'STOP', 'isf_cav: vel_stab gamma formulation had too many iterations ...' )
145            ELSE IF ( zerr <= 0.01_wp ) THEN  ! convergence is achieve
146               lit = .FALSE.
147            ELSE                              ! converge is not yet achieve
148               nit = nit + 1
149               zqh_b(:,:) = zqhc(:,:)+zqoce(:,:)
150            END IF
151         END SELECT
152         !
153      END DO
154      !
155      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
156         ! compute heat and water flux ( > 0 from isf to oce)
157         pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_cav(ji,jj)
158         zqoce(ji,jj) = zqoce(ji,jj) * mskisf_cav(ji,jj)
159         zqhc (ji,jj) = zqhc(ji,jj)  * mskisf_cav(ji,jj)
160         !
161         ! compute heat content flux ( > 0 from isf to oce)
162         zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf    ! 2d latent heat flux (W/m2)
163         !
164         ! total heat flux ( > 0 from isf to oce)
165         zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
166         !
167         ! set temperature content
168         ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
169      END_2D
170      !
171      ! output fluxes
172      CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
173      !
174      ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
175      IF (lrst_oce) CALL isfrst_write(kt, 'cav', ptsc, pqfwf)
176      !
177      IF ( ln_isfdebug ) THEN
178         IF(lwp) WRITE(numout,*) ''
179         CALL debug('isf_cav: ptsc T',ptsc(:,:,1))
180         CALL debug('isf_cav: ptsc S',ptsc(:,:,2))
181         CALL debug('isf_cav: pqfwf fwf',pqfwf(:,:))
182         IF(lwp) WRITE(numout,*) ''
183      END IF
184      !
185   END SUBROUTINE isf_cav
186
187   SUBROUTINE isf_cav_init
188      !!---------------------------------------------------------------------
189      !!                  ***  ROUTINE isf_cav_init ***
190      !!
191      !! ** Purpose : initialisation of variable needed to compute melt under an ice shelf
192      !!
193      !!----------------------------------------------------------------------
194      INTEGER :: ierr
195      !!---------------------------------------------------------------------
196      !
197      !==============
198      ! 0: allocation
199      !==============
200      !
201      CALL isf_alloc_cav()
202      !
203      !==================
204      ! 1: initialisation
205      !==================
206      !
207      ! top and bottom level of the 'top boundary layer'
208      misfkt_cav(:,:)    = mikt(:,:) ; misfkb_cav(:,:)    = 1
209      !
210      ! thickness of 'tbl' and fraction of bottom cell affected by 'tbl'
211      rhisf_tbl_cav(:,:) = 0.0_wp    ; rfrac_tbl_cav(:,:) = 0.0_wp
212      !
213      ! cavity mask
214      mskisf_cav(:,:)    = (1._wp - tmask(:,:,1)) * ssmask(:,:)
215      !================
216      ! 2: activate restart
217      !================
218      !
219      !================
220      ! 3: read restart
221      !================
222      !
223      ! read cav variable from restart
224      IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
225      !
226      !==========================================
227      ! 3: specific allocation and initialisation (depending of scheme choice)
228      !==========================================
229      !
230      SELECT CASE ( TRIM(cn_isfcav_mlt) )
231      CASE( 'spe' )
232
233         ALLOCATE( sf_isfcav_fwf(1), STAT=ierr )
234         ALLOCATE( sf_isfcav_fwf(1)%fnow(jpi,jpj,1), sf_isfcav_fwf(1)%fdta(jpi,jpj,1,2) )
235         CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_isfdir, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
236
237         IF(lwp) WRITE(numout,*)
238         IF(lwp) WRITE(numout,*) '  ==>> The ice shelf melt inside the cavity is read from forcing files'
239
240      CASE( '2eq' )
241         IF(lwp) WRITE(numout,*)
242         IF(lwp) WRITE(numout,*) '  ==>> The original ISOMIP melt formulation is used to compute melt under the ice shelves'
243
244      CASE( '3eq' )
245         ! coeficient for linearisation of potential tfreez
246         ! Crude approximation for pressure (but commonly used)
247         IF ( ln_teos10 ) THEN   ! linearisation from Jourdain et al. (2017)
248            risf_lamb1 =-0.0564_wp
249            risf_lamb2 = 0.0773_wp
250            risf_lamb3 =-7.8633e-8 * grav * rho0
251         ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
252            risf_lamb1 =-0.0573_wp
253            risf_lamb2 = 0.0832_wp
254            risf_lamb3 =-7.5300e-8 * grav * rho0
255         ENDIF
256
257         IF(lwp) WRITE(numout,*)
258         IF(lwp) WRITE(numout,*) '  ==>> The 3 equations melt formulation is used to compute melt under the ice shelves'
259
260      CASE DEFAULT
261         CALL ctl_stop(' cn_isfcav_mlt method unknown (spe, 2eq, 3eq), check namelist')
262      END SELECT
263      !
264   END SUBROUTINE isf_cav_init
265
266END MODULE isfcav
Note: See TracBrowser for help on using the repository browser.