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

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 9406

Last change on this file since 9406 was 9406, checked in by mathiot, 6 years ago
fix for ticket #2066 in nemo 3.6
Property svn:keywords set to `Id`
File size: 45.1 KB

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

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

Download in other formats: