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 @ 5828

Last change on this file since 5828 was 5828, checked in by mathiot, 9 years ago
Fixes for initialisation issues + bug highlited on the MetO machine
Property svn:keywords set to `Id`
File size: 42.0 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 ; fwfisf_b(:,:) = 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	zsfrz=(-zbqe-SQRT(zdis))*zaqer
486	IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
487
488	! compute t freeze
489	zfrz(ji,jj)=zeps4+zlamb1*zsfrz
490
491	! zfwflx is upward water flux
492	! zhtflx is upward heat flux (out of ocean)
493	! compute the upward water and heat flux
494	zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
495	zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (zti(ji,jj) - zfrz(ji,jj) )
496	END DO
497	END DO
498
499	! compute heat and water flux
500	qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
501	fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
502
503	ENDIF
504
505	! define if we need to iterate
506	IF ( nn_gammablk .LT. 2 ) THEN ; lit = .FALSE.
507	ELSE
508	! check total number of iteration
509	IF (nit .GE. 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
510	ELSE ; nit = nit + 1
511	ENDIF
512
513	! compute error between 2 iterations
514	! if needed save gammat and compute zhtflx_b for next iteration
515	zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
516	IF ( zerr .LE. 0.01 ) THEN ; lit = .FALSE.
517	ELSE ; zhtflx_b(:,:) = zhtflx(:,:)
518	ENDIF
519	END IF
520	END DO
521	!
522	! output
523	IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat)
524	IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas)
525	!
526	CALL wrk_dealloc( jpi,jpj, zfrz , zpress, zti , zgammat, zgammas )
527	CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b )
528	!
529	IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav')
530
531	END SUBROUTINE sbc_isf_cav
532
533	SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf )
534	!!----------------------------------------------------------------------
535	!! ** Purpose : compute the coefficient echange for heat flux
536	!!
537	!! ** Method : gamma assume constant or depends of u* and stability
538	!!
539	!! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
540	!! Jenkins et al., 2010, JPO, p2298-2312
541	!!---------------------------------------------------------------------
542	REAL(wp), DIMENSION(:,:), INTENT(out) :: gt, gs
543	REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf
544
545	INTEGER :: ikt ! loop index
546	INTEGER :: ji,jj
547	REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity
548	REAL(wp) :: zdku, zdkv ! U, V shear
549	REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number
550	REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point
551	REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness
552	REAL(wp) :: zhmax ! limitation of mol
553	REAL(wp) :: zetastar ! stability parameter
554	REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence
555	REAL(wp) :: zcoef ! temporary coef
556	REAL(wp) :: zdep
557	REAL(wp), PARAMETER :: zxsiN = 0.052 ! dimensionless constant
558	REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
559	REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
560	REAL(wp) :: zeps = 1.0e-20_wp ! conversion: mm/s ==> m/s
561	REAL(wp), DIMENSION(2) :: zts, zab
562	!!---------------------------------------------------------------------
563	CALL wrk_alloc( jpi,jpj, zustar )
564	!
565	IF ( nn_gammablk == 0 ) THEN
566	!! gamma is constant (specified in namelist)
567	gt(:,:) = rn_gammat0
568	gs(:,:) = rn_gammas0
569
570	ELSE IF ( nn_gammablk == 1 ) THEN
571	!! gamma is assume to be proportional to u*
572	!! Jenkins et al., 2010, JPO, p2298-2312
573
574	!! compute ustar
575	zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
576
577	!! Compute gammats
578	gt(:,:) = zustar(:,:) * rn_gammat0
579	gs(:,:) = zustar(:,:) * rn_gammas0
580
581	ELSE IF ( nn_gammablk == 2 ) THEN
582	!! gamma depends of stability of boundary layer
583	!! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
584	!! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
585	!! compute ustar
586	zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
587
588	!! compute Pr and Sc number (can be improved)
589	zPr = 13.8
590	zSc = 2432.0
591
592	!! compute gamma mole
593	zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
594	zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
595
596	!! compute gamma
597	DO ji=2,jpi
598	DO jj=2,jpj
599	ikt = mikt(ji,jj)
600
601	IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think
602	gt = rn_gammat0
603	gs = rn_gammas0
604	ELSE
605	!! compute Rc number (as done in zdfric.F90)
606	zcoef = 0.5 / fse3w(ji,jj,ikt)
607	! ! shear of horizontal velocity
608	zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) &
609	& -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) )
610	zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) &
611	& -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) )
612	! ! richardson number (minimum value set to zero)
613	zRc = rn2(ji,jj,ikt+1) / MAX( zdkuzdku + zdkvzdkv, zeps )
614
615	!! compute bouyancy
616	zts(jp_tem) = ttbl(ji,jj)
617	zts(jp_sal) = stbl(ji,jj)
618	zdep = fsdepw(ji,jj,ikt)
619	!
620	CALL eos_rab( zts, zdep, zab )
621	!
622	!! compute length scale
623	zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!
624
625	!! compute Monin Obukov Length
626	! Maximum boundary layer depth
627	zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp
628	! Compute Monin obukhov length scale at the surface and Ekman depth:
629	zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + epsln))
630	zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
631
632	!! compute eta* (stability parameter)
633	zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
634
635	!! compute the sublayer thickness
636	zhnu = 5 * znu / zustar(ji,jj)
637
638	!! compute gamma turb
639	zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
640	& + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
641
642	!! compute gammats
643	gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
644	gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
645	END IF
646	END DO
647	END DO
648	CALL lbc_lnk(gt(:,:),'T',1.)
649	CALL lbc_lnk(gs(:,:),'T',1.)
650	END IF
651	CALL wrk_dealloc( jpi,jpj, zustar )
652
653	END SUBROUTINE
654
655	SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
656	!!----------------------------------------------------------------------
657	!! * SUBROUTINE sbc_isf_tbl *
658	!!
659	!! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
660	!!
661	!!----------------------------------------------------------------------
662	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
663	REAL(wp), DIMENSION(:,:) , INTENT(out):: varout
664
665	CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
666
667	REAL(wp) :: ze3, zhk
668	REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl
669
670	INTEGER :: ji,jj,jk
671	INTEGER :: ikt,ikb
672
673	CALL wrk_alloc( jpi,jpj, zhisf_tbl)
674
675	! get first and last level of tbl
676
677	varout(:,:)=0._wp
678	IF (cptin == 'U') THEN
679	DO jj = 1,jpj
680	DO ji = 1,jpi
681	ikt = miku(ji,jj) ; ikb = miku(ji,jj)
682	! thickness of boundary layer at least the top level thickness
683	zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
684
685	! determine the deepest level influenced by the boundary layer
686	! test on tmask useless ?????
687	DO jk = ikt, mbku(ji,jj)
688	IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
689	END DO
690	zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
691
692	! level fully include in the ice shelf boundary layer
693	DO jk = ikt, ikb - 1
694	ze3 = fse3u_n(ji,jj,jk)
695	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
696	END DO
697
698	! level partially include in ice shelf boundary layer
699	zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
700	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
701	END DO
702	END DO
703	DO jj = 2,jpj
704	DO ji = 2,jpi
705	varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj))
706	END DO
707	END DO
708	CALL lbc_lnk(varout,'T',-1.)
709	END IF
710
711	IF (cptin == 'V') THEN
712	DO jj = 1,jpj
713	DO ji = 1,jpi
714	ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
715	! thickness of boundary layer at least the top level thickness
716	zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
717
718	! determine the deepest level influenced by the boundary layer
719	! test on tmask useless ?????
720	DO jk = ikt, mbkv(ji,jj)
721	IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
722	END DO
723	zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
724
725	! level fully include in the ice shelf boundary layer
726	DO jk = ikt, ikb - 1
727	ze3 = fse3v_n(ji,jj,jk)
728	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
729	END DO
730
731	! level partially include in ice shelf boundary layer
732	zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
733	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
734	END DO
735	END DO
736	DO jj = 2,jpj
737	DO ji = 2,jpi
738	varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1))
739	END DO
740	END DO
741	CALL lbc_lnk(varout,'T',-1.)
742	END IF
743
744	IF (cptin == 'T') THEN
745	DO jj = 1,jpj
746	DO ji = 1,jpi
747	ikt = misfkt(ji,jj)
748	ikb = misfkb(ji,jj)
749
750	! level fully include in the ice shelf boundary layer
751	DO jk = ikt, ikb - 1
752	ze3 = fse3t_n(ji,jj,jk)
753	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
754	END DO
755
756	! level partially include in ice shelf boundary layer
757	zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
758	varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
759	END DO
760	END DO
761	END IF
762	varout(:,:) = varout(:,:) * ssmask(:,:)
763
764	CALL wrk_dealloc( jpi,jpj, zhisf_tbl )
765
766	END SUBROUTINE sbc_isf_tbl
767
768
769	SUBROUTINE sbc_isf_div( phdivn )
770	!!----------------------------------------------------------------------
771	!! * SUBROUTINE sbc_isf_div *
772	!!
773	!! ** Purpose : update the horizontal divergence with the runoff inflow
774	!!
775	!! ** Method :
776	!! CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
777	!! divergence and expressed in m/s
778	!!
779	!! ** Action : phdivn decreased by the runoff inflow
780	!!----------------------------------------------------------------------
781	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
782	!!
783	INTEGER :: ji, jj, jk ! dummy loop indices
784	INTEGER :: ikt, ikb
785	REAL(wp) :: zhk
786	REAL(wp) :: zfact ! local scalar
787	!!----------------------------------------------------------------------
788	!
789	zfact = 0.5_wp
790	!
791	IF (lk_vvl) THEN ! need to re compute level distribution of isf fresh water
792	DO jj = 1,jpj
793	DO ji = 1,jpi
794	ikt = misfkt(ji,jj)
795	ikb = misfkt(ji,jj)
796	! thickness of boundary layer at least the top level thickness
797	rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
798
799	! determine the deepest level influenced by the boundary layer
800	! test on tmask useless ?????
801	DO jk = ikt, mbkt(ji,jj)
802	IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
803	END DO
804	rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
805	misfkb(ji,jj) = ikb ! last wet level of the tbl
806	r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
807
808	zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
809	ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer
810	END DO
811	END DO
812	END IF ! vvl case
813	!
814	DO jj = 1,jpj
815	DO ji = 1,jpi
816	ikt = misfkt(ji,jj)
817	ikb = misfkb(ji,jj)
818	! level fully include in the ice shelf boundary layer
819	DO jk = ikt, ikb - 1
820	phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
821	& * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
822	END DO
823	! level partially include in ice shelf boundary layer
824	phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
825	& + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)
826	!== ice shelf melting mass distributed over several levels ==!
827	END DO
828	END DO
829	!
830	END SUBROUTINE sbc_isf_div
831
832	FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
833	!!----------------------------------------------------------------------
834	!! * ROUTINE eos_init *
835	!!
836	!! ** Purpose : Compute the in-situ temperature [Celcius]
837	!!
838	!! ** Method :
839	!!
840	!! Reference : Bryden,h.,1973,deep-sea res.,20,401-408
841	!!----------------------------------------------------------------------
842	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptem ! potential temperature [Celcius]
843	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
844	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ppress ! pressure [dBar]
845	REAL(wp), DIMENSION(:,:), POINTER :: pti ! in-situ temperature [Celcius]
846	REAL(wp) :: zt, zs, zp, zh, zq, zxk
847	INTEGER :: ji, jj ! dummy loop indices
848	!
849	CALL wrk_alloc( jpi,jpj, pti )
850	!
851	DO jj=1,jpj
852	DO ji=1,jpi
853	zh = ppress(ji,jj)
854	! Theta1
855	zt = ptem(ji,jj)
856	zs = psal(ji,jj)
857	zp = 0.0
858	zxk= zh * fsatg( zs, zt, zp )
859	zt = zt + 0.5 * zxk
860	zq = zxk
861	! Theta2
862	zp = zp + 0.5 * zh
863	zxk= zh*fsatg( zs, zt, zp )
864	zt = zt + 0.29289322 * ( zxk - zq )
865	zq = 0.58578644 * zxk + 0.121320344 * zq
866	! Theta3
867	zxk= zh * fsatg( zs, zt, zp )
868	zt = zt + 1.707106781 * ( zxk - zq )
869	zq = 3.414213562 * zxk - 4.121320344 * zq
870	! Theta4
871	zp = zp + 0.5 * zh
872	zxk= zh * fsatg( zs, zt, zp )
873	pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
874	END DO
875	END DO
876	!
877	CALL wrk_dealloc( jpi,jpj, pti )
878	!
879	END FUNCTION tinsitu
880	!
881	FUNCTION fsatg( pfps, pfpt, pfphp )
882	!!----------------------------------------------------------------------
883	!! * FUNCTION fsatg *
884	!!
885	!! ** Purpose : Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
886	!!
887	!! ** Reference : Bryden,h.,1973,deep-sea res.,20,401-408
888	!!
889	!! ** units : pressure pfphp decibars
890	!! temperature pfpt deg celsius (ipts-68)
891	!! salinity pfps (ipss-78)
892	!! adiabatic fsatg deg. c/decibar
893	!!----------------------------------------------------------------------
894	REAL(wp) :: pfps, pfpt, pfphp
895	REAL(wp) :: fsatg
896	!
897	fsatg = (((-2.1687e-16pfpt+1.8676e-14)pfpt-4.6206e-13)*pfphp &
898	& +((2.7759e-12pfpt-1.1351e-10)(pfps-35.)+((-5.4481e-14*pfpt &
899	& +8.733e-12)pfpt-6.7795e-10)pfpt+1.8741e-8))*pfphp &
900	& +(-4.2393e-8pfpt+1.8932e-6)(pfps-35.) &
901	& +((6.6228e-10pfpt-6.836e-8)pfpt+8.5258e-6)*pfpt+3.5803e-5
902	!
903	END FUNCTION fsatg
904	!!======================================================================
905	END MODULE sbcisf

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

Download in other formats: