New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcisf.F90 in branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 8 years ago

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