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/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5075

Last change on this file since 5075 was 5075, checked in by timgraham, 9 years ago

Upgraded branch to current head of trunk (r5072) so it can be used with the trunk

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