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

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

changes in icesheet/iceberg (coupled configuration)

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