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 branches/UKMO/dev_r5518_GO6_package_isfprint/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_isfprint/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 13073

Last change on this file since 13073 was 13073, checked in by timgraham, 5 years ago

Add extra test to some more write statements

File size: 49.6 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 lbclnk          !
21   USE iom             ! I/O manager library
22   USE in_out_manager  ! I/O manager
23   USE wrk_nemo        ! Memory allocation
24   USE timing          ! Timing
25   USE lib_fortran     ! glob_sum
26   USE zdfbfr
27   USE fldread         ! read input field at current time step
28   USE lib_fortran, ONLY: glob_sum
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
34
35   ! public in order to be able to output then
36
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf
39   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
40   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence
41   INTEGER , PUBLIC ::   nn_isfblk                   !:
42   INTEGER , PUBLIC ::   nn_gammablk                 !:
43   LOGICAL , PUBLIC ::   ln_conserve                 !:
44   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
45   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient
46   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
47
48   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
49   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl
50   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl
51   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl
52   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
53   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
54   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
55
56
57   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
58   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
59   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
61   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
62
63!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
64   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
65   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
67   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
69   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
70   
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73   !!----------------------------------------------------------------------
74   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
75   !! $Id$
76   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
77   !!----------------------------------------------------------------------
78
79CONTAINS
80 
81  SUBROUTINE sbc_isf(kt)
82    INTEGER, INTENT(in)          ::   kt         ! ocean time step
83    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
84    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
85    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
86    REAL(wp)                     ::   rmin
87    REAL(wp)                     ::   zhk
88    REAL(wp)                     ::   zt_frz, zpress
89    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
90    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
91    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
92    INTEGER           ::   ios           ! Local integer output status for namelist read
93    LOGICAL                      :: ll_wrtstp  !write diagnostics?
94
95    REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
96    REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
97      !
98      !!---------------------------------------------------------------------
99      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
100                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
101      !
102      !
103      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
104      !                                         ! ====================== !
105      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
106         !                                      ! ====================== !
107         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
108         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
109901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
110
111         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
112         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
113902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
114         IF(lwm .AND. nprint > 2) WRITE ( numond, namsbc_isf )
115
116
117         IF ( lwp ) WRITE(numout,*)
118         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
119         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
120         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
121         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
122         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
123         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
124         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
125         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
126         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
127         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
128         IF(lwp .AND. lflush) CALL flush(numout)
129         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
130            rdivisf = 1._wp
131         ELSE
132            rdivisf = 0._wp
133         END IF
134         !
135         ! Allocate public variable
136         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
137         !
138         ! initialisation
139         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
140         risf_tsc(:,:,:)  = 0._wp
141         !
142         ! define isf tbl tickness, top and bottom indice
143         IF      (nn_isf == 1) THEN
144            rhisf_tbl(:,:) = rn_hisf_tbl
145            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
146         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
147            ALLOCATE( sf_rnfisf(1), STAT=ierror )
148            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
149            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
150
151            !: read effective lenght (BG03)
152            IF (nn_isf == 2) THEN
153               ! Read Data and save some integral values
154               CALL iom_open( sn_Leff_isf%clname, inum )
155               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
156               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
157               CALL iom_close(inum)
158               !
159               risfLeff = risfLeff*1000           !: convertion in m
160            END IF
161
162           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
163            CALL iom_open( sn_depmax_isf%clname, inum )
164            cvarhisf = TRIM(sn_depmax_isf%clvar)
165            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
166            CALL iom_close(inum)
167            !
168            CALL iom_open( sn_depmin_isf%clname, inum )
169            cvarzisf = TRIM(sn_depmin_isf%clvar)
170            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
171            CALL iom_close(inum)
172            !
173            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
174
175           !! compute first level of the top boundary layer
176           DO ji = 1, jpi
177              DO jj = 1, jpj
178                  jk = 2
179                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
180                  misfkt(ji,jj) = jk-1
181               END DO
182            END DO
183
184         ELSE IF ( nn_isf == 4 ) THEN
185            ! as in nn_isf == 1
186            rhisf_tbl(:,:) = rn_hisf_tbl
187            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
188           
189            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
190            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
191            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
192            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
193            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
194            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
195         END IF
196         
197         ! save initial top boundary layer thickness         
198         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
199
200      END IF
201
202      !                                            ! ---------------------------------------- !
203      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
204         !                                         ! ---------------------------------------- !
205         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
206         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
207         !
208      ENDIF
209
210      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
211
212         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
213         DO jj = 1,jpj
214            DO ji = 1,jpi
215               ikt = misfkt(ji,jj)
216               ikb = misfkt(ji,jj)
217               ! thickness of boundary layer at least the top level thickness
218               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
219
220               ! determine the deepest level influenced by the boundary layer
221               DO jk = ikt, mbkt(ji,jj)
222                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
223               END DO
224               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
225               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
226               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
227
228               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
229               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
230            END DO
231         END DO
232
233         ! compute salf and heat flux
234         IF (nn_isf == 1) THEN
235            ! realistic ice shelf formulation
236            ! compute T/S/U/V for the top boundary layer
237            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
238            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
239            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
240            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
241            ! iom print
242            CALL iom_put('ttbl',ttbl(:,:))
243            CALL iom_put('stbl',stbl(:,:))
244            CALL iom_put('utbl',utbl(:,:))
245            CALL iom_put('vtbl',vtbl(:,:))
246            ! compute fwf and heat flux
247            CALL sbc_isf_cav (kt)
248
249         ELSE IF (nn_isf == 2) THEN
250            ! Beckmann and Goosse parametrisation
251            stbl(:,:)   = soce
252            CALL sbc_isf_bg03(kt)
253
254         ELSE IF (nn_isf == 3) THEN
255            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
256            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
257            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
258
259            IF( lk_oasis) THEN
260            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
261            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
262
263              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
264              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
265              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
266 
267              ! All related global sums must be done bit reproducibly
268               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
269
270               ! use ABS function because we need to preserve the sign of fwfisf
271               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
272              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
273              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
274
275               ! check
276               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
277
278               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
279
280               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
281
282               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
283
284               ! use ABS function because we need to preserve the sign of fwfisf
285               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
286              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
287              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
288     
289               ! check
290               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
291
292               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
293
294               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
295
296            ENDIF
297            ENDIF
298
299            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
300            stbl(:,:)   = soce
301
302         ELSE IF (nn_isf == 4) THEN
303            ! specified fwf and heat flux forcing beneath the ice shelf
304            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
305            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
306            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
307
308            IF( lk_oasis) THEN
309            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
310            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
311
312              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
313              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
314              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
315
316              ! All related global sums must be done bit reproducibly
317               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
318
319               ! use ABS function because we need to preserve the sign of fwfisf
320               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
321              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
322              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
323
324               ! check
325               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
326
327               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
328
329               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
330
331               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
332
333               ! use ABS function because we need to preserve the sign of fwfisf
334               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
335              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
336              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
337     
338               ! check
339               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
340
341               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
342
343               IF(lwp .AND. ll_wrtstp .AND. sn_cfctl%l_sbcisf) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
344
345            ENDIF
346            ENDIF
347
348            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
349            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
350            stbl(:,:)   = soce
351
352         END IF
353         ! compute tsc due to isf
354         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
355!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
356         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
357         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
358         
359         ! salt effect already take into account in vertical advection
360         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
361
362         ! output
363         IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf)
364         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce )
365
366         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
367         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
368 
369         ! lbclnk
370         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
371         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
372         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
373         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
374
375!=============================================================================================================================================
376         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
377            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
378            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
379
380            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
381            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
382            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
383            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
384
385            DO jj = 1,jpj
386               DO ji = 1,jpi
387                  ikt = misfkt(ji,jj)
388                  ikb = misfkb(ji,jj)
389                  DO jk = ikt, ikb - 1
390                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
391                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
392                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
393                  END DO
394                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
395                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
396                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
397               END DO
398            END DO
399
400            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
401            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
402            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
403            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
404
405            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
406            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
407         END IF
408!=============================================================================================================================================
409
410         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
411            IF( ln_rstart .AND.    &                     ! Restart: read in restart file
412                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
413               IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
414               IF(lwp .AND. lflush) CALL flush(numout)
415               CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) )   ! before salt content isf_tsc trend
416               CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
417               CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
418            ELSE
419               fwfisf_b(:,:)    = fwfisf(:,:)
420               risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
421            END IF
422         ENDIF
423         !
424      END IF
425 
426  END SUBROUTINE sbc_isf
427
428  INTEGER FUNCTION sbc_isf_alloc()
429      !!----------------------------------------------------------------------
430      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
431      !!----------------------------------------------------------------------
432      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
433      IF( .NOT. ALLOCATED( qisf ) ) THEN
434         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
435               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
436               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
437               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
438               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
439               &    STAT= sbc_isf_alloc )
440         !
441         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
442         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
443         !
444      ENDIF
445  END FUNCTION
446
447  SUBROUTINE sbc_isf_bg03(kt)
448   !!==========================================================================
449   !!                 *** SUBROUTINE sbcisf_bg03  ***
450   !! add net heat and fresh water flux from ice shelf melting
451   !! into the adjacent ocean using the parameterisation by
452   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
453   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
454   !!  (hereafter BG)
455   !!==========================================================================
456   !!----------------------------------------------------------------------
457   !!   sbc_isf_bg03      : routine called from sbcmod
458   !!----------------------------------------------------------------------
459   !!
460   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
461   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
462   !!
463   !! History :
464   !!      !  06-02  (C. Wang) Original code
465   !!----------------------------------------------------------------------
466
467    INTEGER, INTENT ( in ) :: kt
468
469    INTEGER :: ji, jj, jk, jish  !temporary integer
470    INTEGER :: ijkmin
471    INTEGER :: ii, ij, ik 
472    INTEGER :: inum
473
474    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
475    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
476    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
477    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
478   
479    !!----------------------------------------------------------------------
480    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
481     !
482
483    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
484    DO ji = 1, jpi
485       DO jj = 1, jpj
486          ik = misfkt(ji,jj)
487          !! Initialize arrays to 0 (each step)
488          zt_sum = 0.e0_wp
489          IF ( ik .GT. 1 ) THEN
490    ! 3. -----------the average temperature between 200m and 600m ---------------------
491             DO jk = misfkt(ji,jj),misfkb(ji,jj)
492             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
493             ! after verif with UNESCO, wrong sign in BG eq. 2
494             ! Calculate freezing temperature
495                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
496                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
497                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
498             ENDDO
499             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
500   
501    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
502          ! For those corresponding to zonal boundary   
503             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
504                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
505             
506             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                 
507             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
508             !add to salinity trend
509          ELSE
510             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
511          END IF
512       ENDDO
513    ENDDO
514    !
515    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
516  END SUBROUTINE sbc_isf_bg03
517
518   SUBROUTINE sbc_isf_cav( kt )
519      !!---------------------------------------------------------------------
520      !!                     ***  ROUTINE sbc_isf_cav  ***
521      !!
522      !! ** Purpose :   handle surface boundary condition under ice shelf
523      !!
524      !! ** Method  : -
525      !!
526      !! ** Action  :   utau, vtau : remain unchanged
527      !!                taum, wndm : remain unchanged
528      !!                qns        : update heat flux below ice shelf
529      !!                emp, emps  : update freshwater flux below ice shelf
530      !!---------------------------------------------------------------------
531      INTEGER, INTENT(in)          ::   kt         ! ocean time step
532      !
533      LOGICAL :: ln_isomip = .true.
534      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
535      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
536      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
537      REAL(wp) ::   zlamb1, zlamb2, zlamb3
538      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
539      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
540      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
541      REAL(wp) ::   zgammat, zgammas
542      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
543      INTEGER  ::   ji, jj     ! dummy loop indices
544      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
545      INTEGER  ::   ierror     ! return error code
546      LOGICAL  ::   lit=.TRUE.
547      INTEGER  ::   nit
548      !!---------------------------------------------------------------------
549      !
550      ! coeficient for linearisation of tfreez
551      zlamb1=-0.0575
552      zlamb2=0.0901
553      zlamb3=-7.61e-04
554      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
555      !
556      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
557
558      zcfac=0.0_wp 
559      IF (ln_conserve)  zcfac=1.0_wp
560      zpress(:,:)=0.0_wp
561      zgammat2d(:,:)=0.0_wp
562      zgammas2d(:,:)=0.0_wp
563      !
564      !
565!CDIR COLLAPSE
566      DO jj = 1, jpj
567         DO ji = 1, jpi
568            ! Crude approximation for pressure (but commonly used)
569            ! 1e-04 to convert from Pa to dBar
570            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
571            !
572         END DO
573      END DO
574
575! Calculate in-situ temperature (ref to surface)
576      zti(:,:)=tinsitu( ttbl, stbl, zpress )
577! Calculate freezing temperature
578      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
579
580     
581      zhtflx=0._wp ; zfwflx=0._wp
582      IF (nn_isfblk == 1) THEN
583         DO jj = 1, jpj
584            DO ji = 1, jpi
585               IF (mikt(ji,jj) > 1 ) THEN
586                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
587                  DO WHILE ( lit )
588! compute gamma
589                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
590! zhtflx is upward heat flux (out of ocean)
591                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
592! zwflx is upward water flux
593                     zfwflx = - zhtflx/lfusisf
594! test convergence and compute gammat
595                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
596
597                     nit = nit + 1
598                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
599
600! save gammat and compute zhtflx_b
601                     zgammat2d(ji,jj)=zgammat
602                     zhtflx_b = zhtflx
603                  END DO
604
605                  qisf(ji,jj) = - zhtflx
606! For genuine ISOMIP protocol this should probably be something like
607                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
608               ELSE
609                  fwfisf(ji,jj) = 0._wp
610                  qisf(ji,jj)   = 0._wp
611               END IF
612            !
613            END DO
614         END DO
615
616      ELSE IF (nn_isfblk == 2 ) THEN
617
618! More complicated 3 equation thermodynamics as in MITgcm
619!CDIR COLLAPSE
620         DO jj = 2, jpj
621            DO ji = 2, jpi
622               IF (mikt(ji,jj) > 1 ) THEN
623                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
624                  DO WHILE ( lit )
625                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
626
627                     zeps1=rcp*rau0*zgammat
628                     zeps2=lfusisf*rau0*zgammas
629                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
630                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
631                     zeps6=zeps4-zti(ji,jj)
632                     zeps7=zeps4-tsurf
633                     zaqe=zlamb1 * (zeps1 + zeps3)
634                     zaqer=0.5/zaqe
635                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
636                     zcqe=zeps2*stbl(ji,jj)
637                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
638! Presumably zdis can never be negative because gammas is very small compared to gammat
639                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
640                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
641                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
642 
643! zfwflx is upward water flux
644                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
645! zhtflx is upward heat flux (out of ocean)
646! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
647                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
648! zwflx is upward water flux
649! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
650                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
651! test convergence and compute gammat
652                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
653
654                     nit = nit + 1
655                     IF (nit .GE. 51) THEN
656                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
657                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
658                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
659                     END IF
660! save gammat and compute zhtflx_b
661                     zgammat2d(ji,jj)=zgammat
662                     zgammas2d(ji,jj)=zgammas
663                     zhtflx_b = zhtflx
664
665                  END DO
666! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
667                  qisf(ji,jj) = - zhtflx 
668! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
669                  fwfisf(ji,jj) = zfwflx 
670               ELSE
671                  fwfisf(ji,jj) = 0._wp
672                  qisf(ji,jj)   = 0._wp
673               ENDIF
674               !
675            END DO
676         END DO
677      ENDIF
678      ! lbclnk
679      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
680      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
681      ! output
682      CALL iom_put('isfgammat', zgammat2d)
683      CALL iom_put('isfgammas', zgammas2d)
684      !
685      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
686      !
687      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
688
689   END SUBROUTINE sbc_isf_cav
690
691   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
692      !!----------------------------------------------------------------------
693      !! ** Purpose    : compute the coefficient echange for heat flux
694      !!
695      !! ** Method     : gamma assume constant or depends of u* and stability
696      !!
697      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
698      !!                Jenkins et al., 2010, JPO, p2298-2312
699      !!---------------------------------------------------------------------
700      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
701      INTEGER , INTENT(in)    :: ji,jj
702      LOGICAL , INTENT(inout) :: lit
703
704      INTEGER  :: ikt                 ! loop index
705      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
706      REAL(wp) :: zdku, zdkv                 ! U, V shear
707      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
708      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
709      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
710      REAL(wp) :: zhmax                      ! limitation of mol
711      REAL(wp) :: zetastar                   ! stability parameter
712      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
713      REAL(wp) :: zcoef                      ! temporary coef
714      REAL(wp) :: zdep
715      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
716      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
717      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
718      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
719      REAL(wp), DIMENSION(2) :: zts, zab
720      !!---------------------------------------------------------------------
721      !
722      IF( nn_gammablk == 0 ) THEN
723      !! gamma is constant (specified in namelist)
724         gt = rn_gammat0
725         gs = rn_gammas0
726         lit = .FALSE.
727      ELSE IF ( nn_gammablk == 1 ) THEN
728      !! gamma is assume to be proportional to u*
729      !! WARNING in case of Losh 2008 tbl parametrization,
730      !! you have to used the mean value of u in the boundary layer)
731      !! not yet coded
732      !! Jenkins et al., 2010, JPO, p2298-2312
733         ikt = mikt(ji,jj)
734      !! Compute U and V at T points
735   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
736   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
737          zut = utbl(ji,jj)
738          zvt = vtbl(ji,jj)
739
740      !! compute ustar
741         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
742      !! Compute mean value over the TBL
743
744      !! Compute gammats
745         gt = zustar * rn_gammat0
746         gs = zustar * rn_gammas0
747         lit = .FALSE.
748      ELSE IF ( nn_gammablk == 2 ) THEN
749      !! gamma depends of stability of boundary layer
750      !! WARNING in case of Losh 2008 tbl parametrization,
751      !! you have to used the mean value of u in the boundary layer)
752      !! not yet coded
753      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
754      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
755               ikt = mikt(ji,jj)
756
757      !! Compute U and V at T points
758               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
759               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
760
761      !! compute ustar
762               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
763               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
764                 gt = rn_gammat0
765                 gs = rn_gammas0
766               ELSE
767      !! compute Rc number (as done in zdfric.F90)
768               zcoef = 0.5 / fse3w(ji,jj,ikt)
769               !                                            ! shear of horizontal velocity
770               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
771                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
772               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
773                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
774               !                                            ! richardson number (minimum value set to zero)
775               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
776
777      !! compute Pr and Sc number (can be improved)
778               zPr =   13.8
779               zSc = 2432.0
780
781      !! compute gamma mole
782               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
783               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
784
785      !! compute bouyancy
786               zts(jp_tem) = ttbl(ji,jj)
787               zts(jp_sal) = stbl(ji,jj)
788               zdep        = fsdepw(ji,jj,ikt)
789               !
790               CALL eos_rab( zts, zdep, zab )
791                  !
792      !! compute length scale
793               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
794
795      !! compute Monin Obukov Length
796               ! Maximum boundary layer depth
797               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
798               ! Compute Monin obukhov length scale at the surface and Ekman depth:
799               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
800               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
801
802      !! compute eta* (stability parameter)
803               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
804
805      !! compute the sublayer thickness
806               zhnu = 5 * znu / zustar
807      !! compute gamma turb
808               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
809               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
810
811      !! compute gammats
812               gt = zustar / (zgturb + zgmolet)
813               gs = zustar / (zgturb + zgmoles)
814               END IF
815      END IF
816
817   END SUBROUTINE
818
819   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
820      !!----------------------------------------------------------------------
821      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
822      !!
823      !! ** Purpose : compute mean T/S/U/V in the boundary layer
824      !!
825      !!----------------------------------------------------------------------
826      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
827      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
828     
829      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
830
831      REAL(wp) :: ze3, zhk
832      REAL(wp), DIMENSION(:,:), POINTER :: zikt
833
834      INTEGER :: ji,jj,jk
835      INTEGER :: ikt,ikb
836      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
837
838      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
839      CALL wrk_alloc( jpi,jpj, zikt )
840
841      ! get first and last level of tbl
842      mkt(:,:) = misfkt(:,:)
843      mkb(:,:) = misfkb(:,:)
844
845      varout(:,:)=0._wp
846      DO jj = 2,jpj
847         DO ji = 2,jpi
848            IF (ssmask(ji,jj) == 1) THEN
849               ikt = mkt(ji,jj)
850               ikb = mkb(ji,jj)
851
852               ! level fully include in the ice shelf boundary layer
853               DO jk = ikt, ikb - 1
854                  ze3 = fse3t_n(ji,jj,jk)
855                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
856                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
857                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
858                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
859                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
860               END DO
861
862               ! level partially include in ice shelf boundary layer
863               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
864               IF (cptin == 'T') &
865                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
866               IF (cptin == 'U') &
867                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
868               IF (cptin == 'V') &
869                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
870            END IF
871         END DO
872      END DO
873
874      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
875      CALL wrk_dealloc( jpi,jpj, zikt ) 
876
877      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
878      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
879
880   END SUBROUTINE sbc_isf_tbl
881     
882
883   SUBROUTINE sbc_isf_div( phdivn )
884      !!----------------------------------------------------------------------
885      !!                  ***  SUBROUTINE sbc_isf_div  ***
886      !!       
887      !! ** Purpose :   update the horizontal divergence with the runoff inflow
888      !!
889      !! ** Method  :   
890      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
891      !!                          divergence and expressed in m/s
892      !!
893      !! ** Action  :   phdivn   decreased by the runoff inflow
894      !!----------------------------------------------------------------------
895      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
896      !!
897      INTEGER  ::   ji, jj, jk   ! dummy loop indices
898      INTEGER  ::   ikt, ikb 
899      INTEGER  ::   nk_isf
900      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
901      REAL(wp)     ::   zfact     ! local scalar
902      !!----------------------------------------------------------------------
903      !
904      zfact   = 0.5_wp
905      !
906      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
907         DO jj = 1,jpj
908            DO ji = 1,jpi
909               ikt = misfkt(ji,jj)
910               ikb = misfkt(ji,jj)
911               ! thickness of boundary layer at least the top level thickness
912               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
913
914               ! determine the deepest level influenced by the boundary layer
915               ! test on tmask useless ?????
916               DO jk = ikt, mbkt(ji,jj)
917                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
918               END DO
919               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
920               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
921               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
922
923               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
924               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
925            END DO
926         END DO
927      END IF ! vvl case
928      !
929      DO jj = 1,jpj
930         DO ji = 1,jpi
931               ikt = misfkt(ji,jj)
932               ikb = misfkb(ji,jj)
933               ! level fully include in the ice shelf boundary layer
934               DO jk = ikt, ikb - 1
935                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
936                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
937               END DO
938               ! level partially include in ice shelf boundary layer
939               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
940                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
941            !==   ice shelf melting mass distributed over several levels   ==!
942         END DO
943      END DO
944      !
945   END SUBROUTINE sbc_isf_div
946                       
947   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
948      !!----------------------------------------------------------------------
949      !!                 ***  ROUTINE eos_init  ***
950      !!
951      !! ** Purpose :   Compute the in-situ temperature [Celcius]
952      !!
953      !! ** Method  :   
954      !!
955      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
956      !!----------------------------------------------------------------------
957      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
958      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
959      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
960      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
961!      REAL(wp) :: fsatg
962!      REAL(wp) :: pfps, pfpt, pfphp
963      REAL(wp) :: zt, zs, zp, zh, zq, zxk
964      INTEGER  :: ji, jj            ! dummy loop indices
965      !
966      CALL wrk_alloc( jpi,jpj, pti  )
967      !
968      DO jj=1,jpj
969         DO ji=1,jpi
970            zh = ppress(ji,jj)
971! Theta1
972            zt = ptem(ji,jj)
973            zs = psal(ji,jj)
974            zp = 0.0
975            zxk= zh * fsatg( zs, zt, zp )
976            zt = zt + 0.5 * zxk
977            zq = zxk
978! Theta2
979            zp = zp + 0.5 * zh
980            zxk= zh*fsatg( zs, zt, zp )
981            zt = zt + 0.29289322 * ( zxk - zq )
982            zq = 0.58578644 * zxk + 0.121320344 * zq
983! Theta3
984            zxk= zh * fsatg( zs, zt, zp )
985            zt = zt + 1.707106781 * ( zxk - zq )
986            zq = 3.414213562 * zxk - 4.121320344 * zq
987! Theta4
988            zp = zp + 0.5 * zh
989            zxk= zh * fsatg( zs, zt, zp )
990            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
991         END DO
992      END DO
993      !
994      CALL wrk_dealloc( jpi,jpj, pti  )
995      !
996   END FUNCTION tinsitu
997   !
998   FUNCTION fsatg( pfps, pfpt, pfphp )
999      !!----------------------------------------------------------------------
1000      !!                 ***  FUNCTION fsatg  ***
1001      !!
1002      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
1003      !!
1004      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
1005      !!
1006      !! ** units      :   pressure        pfphp    decibars
1007      !!                   temperature     pfpt     deg celsius (ipts-68)
1008      !!                   salinity        pfps     (ipss-78)
1009      !!                   adiabatic       fsatg    deg. c/decibar
1010      !!----------------------------------------------------------------------
1011      REAL(wp) :: pfps, pfpt, pfphp 
1012      REAL(wp) :: fsatg
1013      !
1014      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
1015        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
1016        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
1017        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
1018        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
1019      !
1020    END FUNCTION fsatg
1021    !!======================================================================
1022END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.