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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcisf.F90 in branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5624

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

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

Download in other formats: