source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 19 months ago

GMED 450 add flush after prints

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