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.
sbcisf.F90 in NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90 @ 12576

Last change on this file since 12576 was 12576, checked in by dancopsey, 2 years ago

Merge in iceberg calving stuff from dev_r5518_coupling_GSI7_GSI8_landice from the start of the branch to revision 6023.

File size: 50.3 KB
Line 
1MODULE sbcisf
2   !!======================================================================
3   !!                       ***  MODULE  sbcisf  ***
4   !! Surface module :  update surface ocean boundary condition under ice
5   !!                   shelf
6   !!======================================================================
7   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
8   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
9   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_isf       : update sbc under ice shelf
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE eosbn2         ! equation of state
19   USE sbc_oce        ! surface boundary condition: ocean fields
20   USE zdfdrg         ! vertical physics: top/bottom drag coef.
21   !
22   USE in_out_manager ! I/O manager
23   USE iom            ! I/O library
24   USE fldread        ! read input field at current time step
25   USE lbclnk         !
26   USE lib_fortran    ! glob_sum
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
32
33   ! public in order to be able to output then
34
35   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
36   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified 
37   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
38   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
39   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
40   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
41
42   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s] 
51
52   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
53
54   REAL(wp), PUBLIC, SAVE ::   rcpisf   = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
55   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
56   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
57   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
58   REAL(wp), PUBLIC, SAVE ::   rLfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
59
60!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
61   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
62   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
64   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
66   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
67   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
68   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
69   
70   !!----------------------------------------------------------------------
71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
72   !! $Id$
73   !! Software governed by the CeCILL license (see ./LICENSE)
74   !!----------------------------------------------------------------------
75CONTAINS
76 
77  SUBROUTINE sbc_isf( kt )
78      !!---------------------------------------------------------------------
79      !!                  ***  ROUTINE sbc_isf  ***
80      !!
81      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf
82      !!              melting and freezing
83      !!
84      !! ** Method  :  4 parameterizations are available according to nn_isf
85      !!               nn_isf = 1 : Realistic ice_shelf formulation
86      !!                        2 : Beckmann & Goose parameterization
87      !!                        3 : Specified runoff in deptht (Mathiot & al. )
88      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt   ! ocean time step
91      !
92      INTEGER ::   ji, jj, jk   ! loop index
93      INTEGER ::   ikt, ikb     ! local integers
94      REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 
95      REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep)
96      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d
97      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d
98      !!---------------------------------------------------------------------
99      !
100      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux
101         !
102         SELECT CASE ( nn_isf )
103         CASE ( 1 )    ! realistic ice shelf formulation
104            ! compute T/S/U/V for the top boundary layer
105            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
106            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
107            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
108            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
109            ! iom print
110            CALL iom_put('ttbl',ttbl(:,:))
111            CALL iom_put('stbl',stbl(:,:))
112            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
113            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
114            ! compute fwf and heat flux
115            ! compute fwf and heat flux
116            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
117            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfusisf  ! heat        flux
118            ENDIF
119            !
120         CASE ( 2 )    ! Beckmann and Goosse parametrisation
121            stbl(:,:)   = soce
122            CALL sbc_isf_bg03(kt)
123            !
124         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
125            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
126            IF( .NOT.l_isfcpl ) THEN
127               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
128               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
129            ENDIF
130
131            IF( lk_oasis) THEN
132              ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
133              IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
134
135                ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
136                ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
137                ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
138
139                ! All related global sums must be done bit reproducibly
140                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
141
142                 ! use ABS function because we need to preserve the sign of fwfisf
143                 WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
144                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
145                &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
146
147                 ! check
148                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
149
150                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
151
152                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
153
154                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
155
156                 ! use ABS function because we need to preserve the sign of fwfisf
157                 WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
158                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
159                &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
160
161                 ! check
162                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
163
164                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
165
166                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
167
168              ENDIF
169            ENDIF
170
171            qisf(:,:)   = fwfisf(:,:) * rLfusisf             ! heat flux
172            stbl(:,:)   = soce
173            !
174         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
175            !          ! specified fwf and heat flux forcing beneath the ice shelf
176            IF( .NOT.l_isfcpl ) THEN
177               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
178               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
179               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
180            ENDIF
181
182            IF( lk_oasis) THEN
183              ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
184              IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
185
186                ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
187                ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
188                ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
189
190                ! All related global sums must be done bit reproducibly
191                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
192
193                 ! use ABS function because we need to preserve the sign of fwfisf
194                 WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
195                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
196                &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
197
198                 ! check
199                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
200
201                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
202
203                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
204
205                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
206
207                 ! use ABS function because we need to preserve the sign of fwfisf
208                 WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
209                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
210                &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
211
212                 ! check
213                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
214
215                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
216
217                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
218
219              ENDIF
220            ENDIF
221
222            qisf(:,:)   = fwfisf(:,:) * rLfusisf               ! heat flux
223            stbl(:,:)   = soce
224            !
225         END SELECT
226
227         ! compute tsc due to isf
228         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
229         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
230         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
231         DO jj = 1,jpj
232            DO ji = 1,jpi
233               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
234            END DO
235         END DO
236         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
237         
238         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
239         risf_tsc(:,:,jp_sal) = 0.0_wp
240
241         ! lbclnk
242         CALL lbc_lnk_multi( 'sbcisf', risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.)
243         ! output
244         IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux
245         IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp )   ! isf sensible+latent heat (W/m2)
246         IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat
247         IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign)
248
249         ! Diagnostics
250         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
251            ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) )
252            ALLOCATE( zqhcisf2d(jpi,jpj) )
253            !
254            zfwfisf3d (:,:,:) = 0._wp                         ! 3d ice shelf melting (kg/m2/s)
255            zqhcisf3d (:,:,:) = 0._wp                         ! 3d heat content flux (W/m2)
256            zqlatisf3d(:,:,:) = 0._wp                         ! 3d ice shelf melting latent heat flux (W/m2)
257            zqhcisf2d (:,:)   = fwfisf(:,:) * zt_frz * rcp    ! 2d heat content flux (W/m2)
258            !
259            DO jj = 1,jpj
260               DO ji = 1,jpi
261                  ikt = misfkt(ji,jj)
262                  ikb = misfkb(ji,jj)
263                  DO jk = ikt, ikb - 1
264                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
265                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
266                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
267                  END DO
268                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   & 
269                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
270                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   & 
271                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
272                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   & 
273                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
274               END DO
275            END DO
276            !
277            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
278            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
279            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
280            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
281            !
282            DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d )
283            DEALLOCATE( zqhcisf2d )
284         ENDIF
285         !
286      ENDIF
287
288      IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
289         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
290            &   iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
291            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
292            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:)         , ldxios = lrxios )   ! before salt content isf_tsc trend
293            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content isf_tsc trend
294            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before salt content isf_tsc trend
295         ELSE
296            fwfisf_b(:,:)    = fwfisf(:,:)
297            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
298         ENDIF
299      ENDIF
300      !
301      IF( lrst_oce ) THEN
302         IF(lwp) WRITE(numout,*)
303         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
304            &                    'at it= ', kt,' date= ', ndastp
305         IF(lwp) WRITE(numout,*) '~~~~'
306         IF( lwxios ) CALL iom_swap(      cwxios_context          )
307         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)         , ldxios = lwxios )
308         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios )
309         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios )
310         IF( lwxios ) CALL iom_swap(      cxios_context          )
311      ENDIF
312      !
313   END SUBROUTINE sbc_isf
314
315
316   INTEGER FUNCTION sbc_isf_alloc()
317      !!----------------------------------------------------------------------
318      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
319      !!----------------------------------------------------------------------
320      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
321      IF( .NOT. ALLOCATED( qisf ) ) THEN
322         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
323               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
324               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
325               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
326               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
327               &    STAT= sbc_isf_alloc )
328         !
329         CALL mpp_sum ( 'sbcisf', sbc_isf_alloc )
330         IF( sbc_isf_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_isf_alloc: failed to allocate arrays.' )
331         !
332      ENDIF
333   END FUNCTION
334
335
336  SUBROUTINE sbc_isf_init
337      !!---------------------------------------------------------------------
338      !!                  ***  ROUTINE sbc_isf_init  ***
339      !!
340      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
341      !!
342      !! ** Method  :  4 parameterizations are available according to nn_isf
343      !!               nn_isf = 1 : Realistic ice_shelf formulation
344      !!                        2 : Beckmann & Goose parameterization
345      !!                        3 : Specified runoff in deptht (Mathiot & al. )
346      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
347      !!----------------------------------------------------------------------
348      INTEGER               :: ji, jj, jk           ! loop index
349      INTEGER               :: ik                   ! current level index
350      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
351      INTEGER               :: inum, ierror
352      INTEGER               :: ios                  ! Local integer output status for namelist read
353      REAL(wp)              :: zhk
354      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
355      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
356      !!----------------------------------------------------------------------
357      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
358                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
359      !!----------------------------------------------------------------------
360
361      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
362      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
363901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_isf in reference namelist' )
364
365      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
366      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
367902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist' )
368      IF(lwm) WRITE ( numond, namsbc_isf )
369
370      IF(lwp) WRITE(numout,*)
371      IF(lwp) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
372      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
373      IF(lwp) WRITE(numout,*) '   Namelist namsbc_isf :'
374      IF(lwp) WRITE(numout,*) '      type ice shelf melting/freezing         nn_isf      = ', nn_isf
375      IF(lwp) WRITE(numout,*) '      bulk formulation (nn_isf=1 only)        nn_isfblk   = ', nn_isfblk
376      IF(lwp) WRITE(numout,*) '      thickness of the top boundary layer     rn_hisf_tbl = ', rn_hisf_tbl
377      IF(lwp) WRITE(numout,*) '      gamma formulation                       nn_gammablk = ', nn_gammablk 
378      IF(lwp) WRITE(numout,*) '      gammat coefficient                      rn_gammat0  = ', rn_gammat0 
379      IF(lwp) WRITE(numout,*) '      gammas coefficient                      rn_gammas0  = ', rn_gammas0 
380      IF(lwp) WRITE(numout,*) '      top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top 
381
382
383                           !  1 = presence of ISF    2 = bg03 parametrisation
384                           !  3 = rnf file for isf   4 = ISF fwf specified
385                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
386      !
387      ! Allocate public variable
388      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
389      !
390      ! initialisation
391      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
392      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
393      !
394      ! define isf tbl tickness, top and bottom indice
395      SELECT CASE ( nn_isf )
396      CASE ( 1 ) 
397         IF(lwp) WRITE(numout,*)
398         IF(lwp) WRITE(numout,*) '      ==>>>   presence of under iceshelf seas (nn_isf = 1)'
399         rhisf_tbl(:,:) = rn_hisf_tbl
400         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
401         !
402      CASE ( 2 , 3 )
403         IF( .NOT.l_isfcpl ) THEN
404             ALLOCATE( sf_rnfisf(1), STAT=ierror )
405             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
406             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
407          ENDIF
408          !  read effective lenght (BG03)
409          IF( nn_isf == 2 ) THEN
410            IF(lwp) WRITE(numout,*)
411            IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (nn_isf = 2)'
412            CALL iom_open( sn_Leff_isf%clname, inum )
413            cvarLeff = TRIM(sn_Leff_isf%clvar)
414            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
415            CALL iom_close(inum)
416            !
417            risfLeff = risfLeff*1000.0_wp           !: convertion in m
418         ELSE
419            IF(lwp) WRITE(numout,*)
420            IF(lwp) WRITE(numout,*) '      ==>>>   rnf file for isf (nn_isf = 3)'
421         ENDIF
422         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
423         CALL iom_open( sn_depmax_isf%clname, inum )
424         cvarhisf = TRIM(sn_depmax_isf%clvar)
425         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
426         CALL iom_close(inum)
427         !
428         CALL iom_open( sn_depmin_isf%clname, inum )
429         cvarzisf = TRIM(sn_depmin_isf%clvar)
430         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
431         CALL iom_close(inum)
432         !
433         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
434
435         !! compute first level of the top boundary layer
436         DO ji = 1, jpi
437            DO jj = 1, jpj
438                ik = 2
439!!gm potential bug: use gdepw_0 not _n
440                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
441                misfkt(ji,jj) = ik-1
442            END DO
443         END DO
444         !
445      CASE ( 4 ) 
446         IF(lwp) WRITE(numout,*)
447         IF(lwp) WRITE(numout,*) '      ==>>>   specified fresh water flux in ISF (nn_isf = 4)'
448         ! as in nn_isf == 1
449         rhisf_tbl(:,:) = rn_hisf_tbl
450         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
451         !
452         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
453         IF( .NOT.l_isfcpl ) THEN
454           ALLOCATE( sf_fwfisf(1), STAT=ierror )
455           ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
456           CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
457         ENDIF
458         !
459      CASE DEFAULT
460         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
461      END SELECT
462         
463      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
464
465      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
466      DO jj = 1,jpj
467         DO ji = 1,jpi
468            ikt = misfkt(ji,jj)
469            ikb = misfkt(ji,jj)
470            ! thickness of boundary layer at least the top level thickness
471            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
472
473            ! determine the deepest level influenced by the boundary layer
474            DO jk = ikt+1, mbkt(ji,jj)
475               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
476            END DO
477            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
478            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
479            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
480
481            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
482            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
483         END DO
484      END DO
485
486      IF( lwxios ) THEN
487          CALL iom_set_rstw_var_active('fwf_isf_b')
488          CALL iom_set_rstw_var_active('isf_hc_b')
489          CALL iom_set_rstw_var_active('isf_sc_b')
490      ENDIF
491
492
493  END SUBROUTINE sbc_isf_init
494
495
496  SUBROUTINE sbc_isf_bg03(kt)
497      !!---------------------------------------------------------------------
498      !!                  ***  ROUTINE sbc_isf_bg03  ***
499      !!
500      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
501      !!          into the adjacent ocean
502      !!
503      !! ** Method  :   See reference
504      !!
505      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
506      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
507      !!         (hereafter BG)
508      !! History :  06-02  (C. Wang) Original code
509      !!----------------------------------------------------------------------
510      INTEGER, INTENT ( in ) :: kt
511      !
512      INTEGER  :: ji, jj, jk ! dummy loop index
513      INTEGER  :: ik         ! current level
514      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
515      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
516      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
517      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
518      !!----------------------------------------------------------------------
519      !
520      DO ji = 1, jpi
521         DO jj = 1, jpj
522            ik = misfkt(ji,jj)
523            !! Initialize arrays to 0 (each step)
524            zt_sum = 0.e0_wp
525            IF ( ik > 1 ) THEN
526               ! 1. -----------the average temperature between 200m and 600m ---------------------
527               DO jk = misfkt(ji,jj),misfkb(ji,jj)
528                  ! Calculate freezing temperature
529                  zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04
530                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
531                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
532               END DO
533               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
534               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
535               ! For those corresponding to zonal boundary   
536               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
537                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
538             
539               fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf          !fresh water flux kg/(m2s)                 
540               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
541               !add to salinity trend
542            ELSE
543               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
544            END IF
545         END DO
546      END DO
547      !
548  END SUBROUTINE sbc_isf_bg03
549
550
551  SUBROUTINE sbc_isf_cav( kt )
552      !!---------------------------------------------------------------------
553      !!                     ***  ROUTINE sbc_isf_cav  ***
554      !!
555      !! ** Purpose :   handle surface boundary condition under ice shelf
556      !!
557      !! ** Method  : -
558      !!
559      !! ** Action  :   utau, vtau : remain unchanged
560      !!                taum, wndm : remain unchanged
561      !!                qns        : update heat flux below ice shelf
562      !!                emp, emps  : update freshwater flux below ice shelf
563      !!---------------------------------------------------------------------
564      INTEGER, INTENT(in) ::   kt   ! ocean time step
565      !
566      INTEGER  ::   ji, jj     ! dummy loop indices
567      INTEGER  ::   nit
568      LOGICAL  ::   lit
569      REAL(wp) ::   zlamb1, zlamb2, zlamb3
570      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
571      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
572      REAL(wp) ::   zeps = 1.e-20_wp       
573      REAL(wp) ::   zerr
574      REAL(wp), DIMENSION(jpi,jpj) ::   zfrz
575      REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas 
576      REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b
577      !!---------------------------------------------------------------------
578      !
579      ! coeficient for linearisation of potential tfreez
580      ! Crude approximation for pressure (but commonly used)
581      IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
582         zlamb1 =-0.0564_wp
583         zlamb2 = 0.0773_wp
584         zlamb3 =-7.8633e-8 * grav * rau0
585      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
586         zlamb1 =-0.0573_wp
587         zlamb2 = 0.0832_wp
588         zlamb3 =-7.53e-8 * grav * rau0
589      ENDIF
590      !
591      ! initialisation
592      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
593      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp   
594      zfwflx (:,:) = 0.0_wp
595
596      ! compute ice shelf melting
597      nit = 1 ; lit = .TRUE.
598      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
599         SELECT CASE ( nn_isfblk )
600         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
601            ! Calculate freezing temperature
602            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
603
604            ! compute gammat every where (2d)
605            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
606           
607            ! compute upward heat flux zhtflx and upward water flux zwflx
608            DO jj = 1, jpj
609               DO ji = 1, jpi
610                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
611                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf
612               END DO
613            END DO
614
615            ! Compute heat flux and upward fresh water flux
616            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
617            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
618
619         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
620            ! compute gammat every where (2d)
621            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
622
623            ! compute upward heat flux zhtflx and upward water flux zwflx
624            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
625            DO jj = 1, jpj
626               DO ji = 1, jpi
627                  ! compute coeficient to solve the 2nd order equation
628                  zeps1 = rcp*rau0*zgammat(ji,jj)
629                  zeps2 = rLfusisf*rau0*zgammas(ji,jj)
630                  zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
631                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
632                  zeps6 = zeps4-ttbl(ji,jj)
633                  zeps7 = zeps4-tsurf
634                  zaqe  = zlamb1 * (zeps1 + zeps3)
635                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
636                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
637                  zcqe  = zeps2*stbl(ji,jj)
638                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
639
640                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
641                  ! compute s freeze
642                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
643                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
644
645                  ! compute t freeze (eq. 22)
646                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
647 
648                  ! zfwflx is upward water flux
649                  ! zhtflx is upward heat flux (out of ocean)
650                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
651                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
652                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
653               END DO
654            END DO
655
656            ! compute heat and water flux
657            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
658            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
659
660         END SELECT
661
662         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
663         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
664         ELSE                           
665            ! check total number of iteration
666            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
667            ELSE                 ; nit = nit + 1
668            END IF
669
670            ! compute error between 2 iterations
671            ! if needed save gammat and compute zhtflx_b for next iteration
672            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
673            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
674            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
675            END IF
676         END IF
677      END DO
678      !
679      CALL iom_put('isfgammat', zgammat)
680      CALL iom_put('isfgammas', zgammas)
681      !
682   END SUBROUTINE sbc_isf_cav
683
684
685   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
686      !!----------------------------------------------------------------------
687      !! ** Purpose    : compute the coefficient echange for heat flux
688      !!
689      !! ** Method     : gamma assume constant or depends of u* and stability
690      !!
691      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
692      !!                Jenkins et al., 2010, JPO, p2298-2312
693      !!---------------------------------------------------------------------
694      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !
695      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !
696      !
697      INTEGER  :: ji, jj                     ! loop index
698      INTEGER  :: ikt                        ! local integer
699      REAL(wp) :: zdku, zdkv                 ! U, V shear
700      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
701      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
702      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
703      REAL(wp) :: zhmax                      ! limitation of mol
704      REAL(wp) :: zetastar                   ! stability parameter
705      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
706      REAL(wp) :: zcoef                      ! temporary coef
707      REAL(wp) :: zdep
708      REAL(wp) :: zeps = 1.0e-20_wp   
709      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
710      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
711      REAL(wp), DIMENSION(2) :: zts, zab
712      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
713      !!---------------------------------------------------------------------
714      !
715      SELECT CASE ( nn_gammablk )
716      CASE ( 0 ) ! gamma is constant (specified in namelist)
717         !! ISOMIP formulation (Hunter et al, 2006)
718         pgt(:,:) = rn_gammat0
719         pgs(:,:) = rn_gammas0
720
721      CASE ( 1 ) ! gamma is assume to be proportional to u*
722         !! Jenkins et al., 2010, JPO, p2298-2312
723         !! Adopted by Asay-Davis et al. (2015)
724         !! compute ustar (eq. 24)
725!!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
726!!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
727!! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co
728!!
729!!     ===>>>>    GM  to be done this chrismas
730!!         
731!!gm end
732         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
733
734         !! Compute gammats
735         pgt(:,:) = zustar(:,:) * rn_gammat0
736         pgs(:,:) = zustar(:,:) * rn_gammas0
737     
738      CASE ( 2 ) ! gamma depends of stability of boundary layer
739         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
740         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
741         !! compute ustar
742         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
743
744         !! compute Pr and Sc number (can be improved)
745         zPr =   13.8_wp
746         zSc = 2432.0_wp
747
748         !! compute gamma mole
749         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
750         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
751
752         !! compute gamma
753         DO ji = 2, jpi
754            DO jj = 2, jpj
755               ikt = mikt(ji,jj)
756
757               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
758                  pgt = rn_gammat0
759                  pgs = rn_gammas0
760               ELSE
761                  !! compute Rc number (as done in zdfric.F90)
762!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
763!!gm moreover, use Max(rn2,0) to take care of static instabilities....
764                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
765                  !                                            ! shear of horizontal velocity
766                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
767                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
768                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
769                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
770                  !                                            ! richardson number (minimum value set to zero)
771                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
772
773                  !! compute bouyancy
774                  zts(jp_tem) = ttbl(ji,jj)
775                  zts(jp_sal) = stbl(ji,jj)
776                  zdep        = gdepw_n(ji,jj,ikt)
777                  !
778                  CALL eos_rab( zts, zdep, zab )
779                  !
780                  !! compute length scale
781                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
782
783                  !! compute Monin Obukov Length
784                  ! Maximum boundary layer depth
785                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
786                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
787                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
788                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
789
790                  !! compute eta* (stability parameter)
791                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
792
793                  !! compute the sublayer thickness
794                  zhnu = 5 * znu / zustar(ji,jj)
795
796                  !! compute gamma turb
797                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
798                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
799
800                  !! compute gammats
801                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
802                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
803               END IF
804            END DO
805         END DO
806         CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.)
807      END SELECT
808      !
809   END SUBROUTINE sbc_isf_gammats
810
811
812   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
813      !!----------------------------------------------------------------------
814      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
815      !!
816      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
817      !!
818      !!----------------------------------------------------------------------
819      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
820      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
821      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
822      !
823      INTEGER ::   ji, jj, jk                ! loop index
824      INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl
825      REAL(wp) ::   ze3, zhk
826      REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
827      !!----------------------------------------------------------------------
828     
829      ! initialisation
830      pvarout(:,:)=0._wp
831   
832      SELECT CASE ( cd_ptin )
833      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
834         DO jj = 1,jpj
835            DO ji = 1,jpi
836               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
837               ! thickness of boundary layer at least the top level thickness
838               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
839
840               ! determine the deepest level influenced by the boundary layer
841               DO jk = ikt+1, mbku(ji,jj)
842                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
843               END DO
844               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
845
846               ! level fully include in the ice shelf boundary layer
847               DO jk = ikt, ikb - 1
848                  ze3 = e3u_n(ji,jj,jk)
849                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
850               END DO
851
852               ! level partially include in ice shelf boundary layer
853               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
854               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
855            END DO
856         END DO
857         DO jj = 2, jpj
858            DO ji = 2, jpi
859!!gm a wet-point only average should be used here !!!
860               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
861            END DO
862         END DO
863         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
864     
865      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
866         DO jj = 1,jpj
867            DO ji = 1,jpi
868               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
869               ! thickness of boundary layer at least the top level thickness
870               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
871
872               ! determine the deepest level influenced by the boundary layer
873               DO jk = ikt+1, mbkv(ji,jj)
874                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
875               END DO
876               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
877
878               ! level fully include in the ice shelf boundary layer
879               DO jk = ikt, ikb - 1
880                  ze3 = e3v_n(ji,jj,jk)
881                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
882               END DO
883
884               ! level partially include in ice shelf boundary layer
885               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
886               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
887            END DO
888         END DO
889         DO jj = 2, jpj
890            DO ji = 2, jpi
891!!gm a wet-point only average should be used here !!!
892               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
893            END DO
894         END DO
895         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
896
897      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
898         DO jj = 1,jpj
899            DO ji = 1,jpi
900               ikt = misfkt(ji,jj)
901               ikb = misfkb(ji,jj)
902
903               ! level fully include in the ice shelf boundary layer
904               DO jk = ikt, ikb - 1
905                  ze3 = e3t_n(ji,jj,jk)
906                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
907               END DO
908
909               ! level partially include in ice shelf boundary layer
910               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
911               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
912            END DO
913         END DO
914      END SELECT
915      !
916      ! mask mean tbl value
917      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
918      !
919   END SUBROUTINE sbc_isf_tbl
920     
921
922   SUBROUTINE sbc_isf_div( phdivn )
923      !!----------------------------------------------------------------------
924      !!                  ***  SUBROUTINE sbc_isf_div  ***
925      !!       
926      !! ** Purpose :   update the horizontal divergence with the runoff inflow
927      !!
928      !! ** Method  :   
929      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
930      !!                          divergence and expressed in m/s
931      !!
932      !! ** Action  :   phdivn   decreased by the runoff inflow
933      !!----------------------------------------------------------------------
934      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
935      !
936      INTEGER  ::   ji, jj, jk   ! dummy loop indices
937      INTEGER  ::   ikt, ikb 
938      REAL(wp) ::   zhk
939      REAL(wp) ::   zfact     ! local scalar
940      !!----------------------------------------------------------------------
941      !
942      zfact   = 0.5_wp
943      !
944      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
945         DO jj = 1,jpj
946            DO ji = 1,jpi
947               ikt = misfkt(ji,jj)
948               ikb = misfkt(ji,jj)
949               ! thickness of boundary layer at least the top level thickness
950               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
951
952               ! determine the deepest level influenced by the boundary layer
953               DO jk = ikt, mbkt(ji,jj)
954                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
955               END DO
956               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
957               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
958               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
959
960               zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
961               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
962            END DO
963         END DO
964      END IF 
965      !
966      !==   ice shelf melting distributed over several levels   ==!
967      DO jj = 1,jpj
968         DO ji = 1,jpi
969               ikt = misfkt(ji,jj)
970               ikb = misfkb(ji,jj)
971               ! level fully include in the ice shelf boundary layer
972               DO jk = ikt, ikb - 1
973                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
974                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
975               END DO
976               ! level partially include in ice shelf boundary layer
977               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
978                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
979         END DO
980      END DO
981      !
982   END SUBROUTINE sbc_isf_div
983
984   !!======================================================================
985END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.