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

Context Navigation

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

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago
GMED 450 add flush after prints
File size: 49.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: