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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5866

Last change on this file since 5866 was 5866, checked in by gm, 8 years ago

#1613: vvl by default: add ln_linssh and remove key_vvl

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