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 NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90 @ 12576

Last change on this file since 12576 was 12576, checked in by dancopsey, 4 years ago
Merge in iceberg calving stuff from dev_r5518_coupling_GSI7_GSI8_landice from the start of the branch to revision 6023.
File size: 50.3 KB

Rev	Line
[4666]	1	MODULE sbcisf
	2	!!======================================================================
	3	!! * MODULE sbcisf *
	4	!! Surface module : update surface ocean boundary condition under ice
	5	!! shelf
	6	!!======================================================================
[9019]	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
[4666]	10	!!----------------------------------------------------------------------
	11
	12	!!----------------------------------------------------------------------
[9019]	13	!! sbc_isf : update sbc under ice shelf
[4666]	14	!!----------------------------------------------------------------------
[9019]	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 zdfdrg ! vertical physics: top/bottom drag coef.
[5836]	21	!
[9019]	22	USE in_out_manager ! I/O manager
	23	USE iom ! I/O library
	24	USE fldread ! read input field at current time step
	25	USE lbclnk !
	26	USE lib_fortran ! glob_sum
[4666]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[7816]	31	PUBLIC sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and divhor
[4666]	32
	33	! public in order to be able to output then
	34
	35	REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m]
[6140]	36	INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified
	37	INTEGER , PUBLIC :: nn_isfblk !: flag to choose the bulk formulation to compute the ice shelf melting
	38	INTEGER , PUBLIC :: nn_gammablk !: flag to choose how the exchange coefficient is computed
	39	REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient []
	40	REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient []
[4666]	41
[9019]	42	INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt , misfkb !: Level of ice shelf base
	43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rzisf_tbl !: depth of calving front (shallowest point) nn_isf ==2/3
	44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl, rhisf_tbl_0 !: thickness of tbl [m]
	45	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hisf_tbl !: 1/thickness of tbl
	46	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ralpha !: proportion of bottom cell influenced by tbl
	47	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !: effective length (Leff) BG03 nn_isf==2
	48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ttbl, stbl, utbl, vtbl !: top boundary layer variable at T point
	49	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2]
	50	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s]
[4666]	51
[9019]	52	LOGICAL, PUBLIC :: l_isfcpl = .false. !: isf recieved from oasis
[7788]	53
[9935]	54	REAL(wp), PUBLIC, SAVE :: rcpisf = 2000.0_wp !: specific heat of ice shelf [J/kg/K]
[9019]	55	REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp !: heat diffusivity through the ice-shelf [m2/s]
	56	REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3]
	57	REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp !: air temperature on top of ice shelf [C]
[9935]	58	REAL(wp), PUBLIC, SAVE :: rLfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg]
[4666]	59
	60	!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
[6140]	61	CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files
	62	TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf melting file to be read
	63	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
	64	TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the isf melting param. file to be read
	65	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf
	66	TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the grounding line depth file to be read
	67	TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the calving line depth file to be read
	68	TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the effective length file to be read
[4666]	69
	70	!!----------------------------------------------------------------------
[9598]	71	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]	72	!! $Id$
[10068]	73	!! Software governed by the CeCILL license (see ./LICENSE)
[4666]	74	!!----------------------------------------------------------------------
	75	CONTAINS
	76
[9019]	77	SUBROUTINE sbc_isf( kt )
[5836]	78	!!---------------------------------------------------------------------
[6140]	79	!! * ROUTINE sbc_isf *
	80	!!
	81	!! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf
	82	!! melting and freezing
	83	!!
	84	!! ** Method : 4 parameterizations are available according to nn_isf
	85	!! nn_isf = 1 : Realistic ice_shelf formulation
	86	!! 2 : Beckmann & Goose parameterization
	87	!! 3 : Specified runoff in deptht (Mathiot & al. )
	88	!! 4 : specified fwf and heat flux forcing beneath the ice shelf
	89	!!----------------------------------------------------------------------
[9019]	90	INTEGER, INTENT(in) :: kt ! ocean time step
[4666]	91	!
[9019]	92	INTEGER :: ji, jj, jk ! loop index
	93	INTEGER :: ikt, ikb ! local integers
[12576]	94	REAL(wp) :: zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
[9019]	95	REAL(wp), DIMENSION(jpi,jpj) :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)
	96	REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zqhcisf2d
	97	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfwfisf3d, zqhcisf3d, zqlatisf3d
[4666]	98	!!---------------------------------------------------------------------
	99	!
[9019]	100	IF( MOD( kt-1, nn_fsbc) == 0 ) THEN ! compute salt and heat flux
	101	!
[6140]	102	SELECT CASE ( nn_isf )
	103	CASE ( 1 ) ! realistic ice shelf formulation
[4666]	104	! compute T/S/U/V for the top boundary layer
	105	CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
	106	CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
[6140]	107	CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')
	108	CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')
[4666]	109	! iom print
	110	CALL iom_put('ttbl',ttbl(:,:))
	111	CALL iom_put('stbl',stbl(:,:))
[6140]	112	CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
	113	CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
[4666]	114	! compute fwf and heat flux
[7788]	115	! compute fwf and heat flux
	116	IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt)
[9935]	117	ELSE ; qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux
[7788]	118	ENDIF
[9019]	119	!
[6140]	120	CASE ( 2 ) ! Beckmann and Goosse parametrisation
[4666]	121	stbl(:,:) = soce
	122	CALL sbc_isf_bg03(kt)
[9019]	123	!
[6140]	124	CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation)
[7788]	125	! specified runoff in depth (Mathiot et al., XXXX in preparation)
	126	IF( .NOT.l_isfcpl ) THEN
	127	CALL fld_read ( kt, nn_fsbc, sf_rnfisf )
	128	fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting)
	129	ENDIF
[12576]	130
	131	IF( lk_oasis) THEN
	132	! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
	133	IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
	134
	135	! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
	136	! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
	137	! to preserve total freshwater conservation in coupled models without an active ice sheet model.
	138
	139	! All related global sums must be done bit reproducibly
	140	zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
	141
	142	! use ABS function because we need to preserve the sign of fwfisf
	143	WHERE( greenland_icesheet_mask(:,:) == 1.0 ) &
	144	& fwfisf(:,:) = fwfisf(:,:) * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
	145	& / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
	146
	147	! check
	148	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
	149
	150	zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
	151
	152	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
	153
	154	zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
	155
	156	! use ABS function because we need to preserve the sign of fwfisf
	157	WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
	158	& fwfisf(:,:) = fwfisf(:,:) * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
	159	& / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
	160
	161	! check
	162	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
	163
	164	zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
	165
	166	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
	167
	168	ENDIF
	169	ENDIF
	170
[9935]	171	qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux
[4666]	172	stbl(:,:) = soce
[9019]	173	!
[6140]	174	CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf
[9019]	175	! ! specified fwf and heat flux forcing beneath the ice shelf
[7788]	176	IF( .NOT.l_isfcpl ) THEN
	177	CALL fld_read ( kt, nn_fsbc, sf_fwfisf )
	178	!CALL fld_read ( kt, nn_fsbc, sf_qisf )
	179	fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1) ! fwf
	180	ENDIF
[12576]	181
	182	IF( lk_oasis) THEN
	183	! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
	184	IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
	185
	186	! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
	187	! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
	188	! to preserve total freshwater conservation in coupled models without an active ice sheet model.
	189
	190	! All related global sums must be done bit reproducibly
	191	zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
	192
	193	! use ABS function because we need to preserve the sign of fwfisf
	194	WHERE( greenland_icesheet_mask(:,:) == 1.0 ) &
	195	& fwfisf(:,:) = fwfisf(:,:) * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
	196	& / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
	197
	198	! check
	199	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
	200
	201	zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
	202
	203	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
	204
	205	zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
	206
	207	! use ABS function because we need to preserve the sign of fwfisf
	208	WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
	209	& fwfisf(:,:) = fwfisf(:,:) * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
	210	& / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
	211
	212	! check
	213	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
	214
	215	zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
	216
	217	IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
	218
	219	ENDIF
	220	ENDIF
	221
[9935]	222	qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux
[4666]	223	stbl(:,:) = soce
[9019]	224	!
[6140]	225	END SELECT
	226
[4666]	227	! compute tsc due to isf
[6140]	228	! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
	229	! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
	230	! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
	231	DO jj = 1,jpj
	232	DO ji = 1,jpi
	233	zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
	234	END DO
	235	END DO
	236	CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
[4666]	237
[6140]	238	risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
	239	risf_tsc(:,:,jp_sal) = 0.0_wp
[5643]	240
[4666]	241	! lbclnk
[10425]	242	CALL lbc_lnk_multi( 'sbcisf', risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.)
[7788]	243	! output
[7968]	244	IF( iom_use('iceshelf_cea') ) CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) ) ! isf mass flux
	245	IF( iom_use('hflx_isf_cea') ) CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp ) ! isf sensible+latent heat (W/m2)
	246	IF( iom_use('qlatisf' ) ) CALL iom_put( 'qlatisf' , qisf(:,:) ) ! isf latent heat
	247	IF( iom_use('fwfisf' ) ) CALL iom_put( 'fwfisf' , fwfisf(:,:) ) ! isf mass flux (opposite sign)
[7788]	248
[9019]	249	! Diagnostics
	250	IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
	251	ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) )
	252	ALLOCATE( zqhcisf2d(jpi,jpj) )
	253	!
	254	zfwfisf3d (:,:,:) = 0._wp ! 3d ice shelf melting (kg/m2/s)
	255	zqhcisf3d (:,:,:) = 0._wp ! 3d heat content flux (W/m2)
	256	zqlatisf3d(:,:,:) = 0._wp ! 3d ice shelf melting latent heat flux (W/m2)
	257	zqhcisf2d (:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2)
	258	!
[7646]	259	DO jj = 1,jpj
	260	DO ji = 1,jpi
	261	ikt = misfkt(ji,jj)
	262	ikb = misfkb(ji,jj)
	263	DO jk = ikt, ikb - 1
	264	zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
	265	zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
	266	zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
	267	END DO
[8329]	268	zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) &
	269	& * ralpha(ji,jj) * e3t_n(ji,jj,jk)
	270	zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) &
	271	& * ralpha(ji,jj) * e3t_n(ji,jj,jk)
	272	zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) &
	273	& * ralpha(ji,jj) * e3t_n(ji,jj,jk)
[7646]	274	END DO
	275	END DO
[9019]	276	!
[7646]	277	CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
	278	CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
	279	CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
	280	CALL iom_put('qhcisf' , zqhcisf2d (:,: ))
[9019]	281	!
	282	DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d )
	283	DEALLOCATE( zqhcisf2d )
[7788]	284	ENDIF
	285	!
[9019]	286	ENDIF
[4666]	287
[9019]	288	IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
	289	IF( ln_rstart .AND. & ! Restart: read in restart file
	290	& iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
	291	IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file'
[9367]	292	CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) , ldxios = lrxios ) ! before salt content isf_tsc trend
	293	CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios ) ! before salt content isf_tsc trend
	294	CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios ) ! before salt content isf_tsc trend
[9019]	295	ELSE
	296	fwfisf_b(:,:) = fwfisf(:,:)
	297	risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
	298	ENDIF
	299	ENDIF
	300	!
	301	IF( lrst_oce ) THEN
	302	IF(lwp) WRITE(numout,*)
	303	IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', &
	304	& 'at it= ', kt,' date= ', ndastp
	305	IF(lwp) WRITE(numout,*) '~~~~'
[9367]	306	IF( lwxios ) CALL iom_swap( cwxios_context )
	307	CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) , ldxios = lwxios )
	308	CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios )
	309	CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios )
	310	IF( lwxios ) CALL iom_swap( cxios_context )
[9019]	311	ENDIF
	312	!
	313	END SUBROUTINE sbc_isf
[5836]	314
[9019]	315
	316	INTEGER FUNCTION sbc_isf_alloc()
[4666]	317	!!----------------------------------------------------------------------
	318	!! * FUNCTION sbc_isf_rnf_alloc *
	319	!!----------------------------------------------------------------------
	320	sbc_isf_alloc = 0 ! set to zero if no array to be allocated
	321	IF( .NOT. ALLOCATED( qisf ) ) THEN
[5120]	322	ALLOCATE( risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj) , &
	323	& rhisf_tbl(jpi,jpj) , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj) , &
	324	& ttbl(jpi,jpj) , stbl(jpi,jpj) , utbl(jpi,jpj) , &
	325	& vtbl(jpi, jpj) , risfLeff(jpi,jpj) , rhisf_tbl_0(jpi,jpj), &
	326	& ralpha(jpi,jpj) , misfkt(jpi,jpj) , misfkb(jpi,jpj) , &
[4946]	327	& STAT= sbc_isf_alloc )
[4666]	328	!
[10425]	329	CALL mpp_sum ( 'sbcisf', sbc_isf_alloc )
	330	IF( sbc_isf_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf_alloc: failed to allocate arrays.' )
[4666]	331	!
[9019]	332	ENDIF
	333	END FUNCTION
[4666]	334
[9019]	335
[6140]	336	SUBROUTINE sbc_isf_init
	337	!!---------------------------------------------------------------------
	338	!! * ROUTINE sbc_isf_init *
	339	!!
	340	!! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
	341	!!
	342	!! ** Method : 4 parameterizations are available according to nn_isf
	343	!! nn_isf = 1 : Realistic ice_shelf formulation
	344	!! 2 : Beckmann & Goose parameterization
	345	!! 3 : Specified runoff in deptht (Mathiot & al. )
	346	!! 4 : specified fwf and heat flux forcing beneath the ice shelf
[5836]	347	!!----------------------------------------------------------------------
[6140]	348	INTEGER :: ji, jj, jk ! loop index
	349	INTEGER :: ik ! current level index
	350	INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer
	351	INTEGER :: inum, ierror
	352	INTEGER :: ios ! Local integer output status for namelist read
	353	REAL(wp) :: zhk
	354	CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file
	355	CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale
[5836]	356	!!----------------------------------------------------------------------
[6140]	357	NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
	358	& sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
	359	!!----------------------------------------------------------------------
	360
	361	REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs
	362	READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
[11536]	363	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist' )
[6140]	364
	365	REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs
	366	READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
[11536]	367	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist' )
[6140]	368	IF(lwm) WRITE ( numond, namsbc_isf )
	369
[9168]	370	IF(lwp) WRITE(numout,*)
	371	IF(lwp) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
	372	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	373	IF(lwp) WRITE(numout,*) ' Namelist namsbc_isf :'
	374	IF(lwp) WRITE(numout,*) ' type ice shelf melting/freezing nn_isf = ', nn_isf
	375	IF(lwp) WRITE(numout,*) ' bulk formulation (nn_isf=1 only) nn_isfblk = ', nn_isfblk
	376	IF(lwp) WRITE(numout,*) ' thickness of the top boundary layer rn_hisf_tbl = ', rn_hisf_tbl
	377	IF(lwp) WRITE(numout,*) ' gamma formulation nn_gammablk = ', nn_gammablk
	378	IF(lwp) WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0
	379	IF(lwp) WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0
	380	IF(lwp) WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top
	381
	382
	383	! 1 = presence of ISF 2 = bg03 parametrisation
	384	! 3 = rnf file for isf 4 = ISF fwf specified
	385	! option 1 and 4 need ln_isfcav = .true. (domzgr)
[6140]	386	!
	387	! Allocate public variable
	388	IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
	389	!
	390	! initialisation
[9019]	391	qisf (:,:) = 0._wp ; fwfisf (:,:) = 0._wp
	392	risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp
[6140]	393	!
	394	! define isf tbl tickness, top and bottom indice
	395	SELECT CASE ( nn_isf )
	396	CASE ( 1 )
[9168]	397	IF(lwp) WRITE(numout,*)
	398	IF(lwp) WRITE(numout,*) ' ==>>> presence of under iceshelf seas (nn_isf = 1)'
[6140]	399	rhisf_tbl(:,:) = rn_hisf_tbl
[9019]	400	misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv
[9168]	401	!
[6140]	402	CASE ( 2 , 3 )
[7788]	403	IF( .NOT.l_isfcpl ) THEN
	404	ALLOCATE( sf_rnfisf(1), STAT=ierror )
	405	ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
	406	CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
	407	ENDIF
	408	! read effective lenght (BG03)
[9168]	409	IF( nn_isf == 2 ) THEN
	410	IF(lwp) WRITE(numout,*)
	411	IF(lwp) WRITE(numout,*) ' ==>>> bg03 parametrisation (nn_isf = 2)'
[6140]	412	CALL iom_open( sn_Leff_isf%clname, inum )
	413	cvarLeff = TRIM(sn_Leff_isf%clvar)
	414	CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
	415	CALL iom_close(inum)
	416	!
	417	risfLeff = risfLeff*1000.0_wp !: convertion in m
[9168]	418	ELSE
	419	IF(lwp) WRITE(numout,*)
	420	IF(lwp) WRITE(numout,*) ' ==>>> rnf file for isf (nn_isf = 3)'
	421	ENDIF
[6140]	422	! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
	423	CALL iom_open( sn_depmax_isf%clname, inum )
	424	cvarhisf = TRIM(sn_depmax_isf%clvar)
	425	CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
	426	CALL iom_close(inum)
	427	!
	428	CALL iom_open( sn_depmin_isf%clname, inum )
	429	cvarzisf = TRIM(sn_depmin_isf%clvar)
	430	CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
	431	CALL iom_close(inum)
	432	!
	433	rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer
	434
	435	!! compute first level of the top boundary layer
	436	DO ji = 1, jpi
	437	DO jj = 1, jpj
	438	ik = 2
[9019]	439	!!gm potential bug: use gdepw_0 not _n
[6140]	440	DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO
	441	misfkt(ji,jj) = ik-1
	442	END DO
	443	END DO
[9168]	444	!
[6140]	445	CASE ( 4 )
[9168]	446	IF(lwp) WRITE(numout,*)
	447	IF(lwp) WRITE(numout,*) ' ==>>> specified fresh water flux in ISF (nn_isf = 4)'
[6140]	448	! as in nn_isf == 1
	449	rhisf_tbl(:,:) = rn_hisf_tbl
[9019]	450	misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv
[9168]	451	!
[6140]	452	! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
[7788]	453	IF( .NOT.l_isfcpl ) THEN
	454	ALLOCATE( sf_fwfisf(1), STAT=ierror )
	455	ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
	456	CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
	457	ENDIF
[9168]	458	!
	459	CASE DEFAULT
	460	CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
[6140]	461	END SELECT
	462
	463	rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
	464
	465	! compute bottom level of isf tbl and thickness of tbl below the ice shelf
	466	DO jj = 1,jpj
	467	DO ji = 1,jpi
	468	ikt = misfkt(ji,jj)
	469	ikb = misfkt(ji,jj)
	470	! thickness of boundary layer at least the top level thickness
	471	rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
	472
	473	! determine the deepest level influenced by the boundary layer
	474	DO jk = ikt+1, mbkt(ji,jj)
[9019]	475	IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
[6140]	476	END DO
	477	rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
	478	misfkb(ji,jj) = ikb ! last wet level of the tbl
	479	r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
	480
	481	zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
	482	ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer
	483	END DO
	484	END DO
	485
[9367]	486	IF( lwxios ) THEN
	487	CALL iom_set_rstw_var_active('fwf_isf_b')
	488	CALL iom_set_rstw_var_active('isf_hc_b')
	489	CALL iom_set_rstw_var_active('isf_sc_b')
	490	ENDIF
	491
	492
[6140]	493	END SUBROUTINE sbc_isf_init
	494
[9019]	495
[6140]	496	SUBROUTINE sbc_isf_bg03(kt)
	497	!!---------------------------------------------------------------------
	498	!! * ROUTINE sbc_isf_bg03 *
[5836]	499	!!
[6140]	500	!! ** Purpose : add net heat and fresh water flux from ice shelf melting
	501	!! into the adjacent ocean
[5836]	502	!!
[6140]	503	!! ** Method : See reference
	504	!!
	505	!! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
	506	!! interaction for climate models", Ocean Modelling 5(2003) 157-170.
	507	!! (hereafter BG)
[9019]	508	!! History : 06-02 (C. Wang) Original code
[5836]	509	!!----------------------------------------------------------------------
	510	INTEGER, INTENT ( in ) :: kt
	511	!
[6140]	512	INTEGER :: ji, jj, jk ! dummy loop index
	513	INTEGER :: ik ! current level
	514	REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m
	515	REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m
	516	REAL(wp) :: zt_frz ! freezing point temperature at depth z
	517	REAL(wp) :: zpress ! pressure to compute the freezing point in depth
	518	!!----------------------------------------------------------------------
	519	!
	520	DO ji = 1, jpi
	521	DO jj = 1, jpj
	522	ik = misfkt(ji,jj)
	523	!! Initialize arrays to 0 (each step)
	524	zt_sum = 0.e0_wp
	525	IF ( ik > 1 ) THEN
	526	! 1. -----------the average temperature between 200m and 600m ---------------------
	527	DO jk = misfkt(ji,jj),misfkb(ji,jj)
	528	! Calculate freezing temperature
[9865]	529	zpress = gravrau0gdept_n(ji,jj,ik)*1.e-04
[6140]	530	CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)
	531	zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp
	532	END DO
	533	zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
	534	! 2. ------------Net heat flux and fresh water flux due to the ice shelf
	535	! For those corresponding to zonal boundary
	536	qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave &
	537	& * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
[4666]	538
[9935]	539	fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf !fresh water flux kg/(m2s)
[6140]	540	fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
	541	!add to salinity trend
	542	ELSE
[9019]	543	qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
[6140]	544	END IF
	545	END DO
	546	END DO
[5836]	547	!
[4666]	548	END SUBROUTINE sbc_isf_bg03
	549
[9019]	550
[6140]	551	SUBROUTINE sbc_isf_cav( kt )
[4666]	552	!!---------------------------------------------------------------------
	553	!! * ROUTINE sbc_isf_cav *
	554	!!
	555	!! ** Purpose : handle surface boundary condition under ice shelf
	556	!!
	557	!! ** Method : -
	558	!!
	559	!! ** Action : utau, vtau : remain unchanged
	560	!! taum, wndm : remain unchanged
	561	!! qns : update heat flux below ice shelf
	562	!! emp, emps : update freshwater flux below ice shelf
	563	!!---------------------------------------------------------------------
[9019]	564	INTEGER, INTENT(in) :: kt ! ocean time step
[4666]	565	!
[6140]	566	INTEGER :: ji, jj ! dummy loop indices
	567	INTEGER :: nit
[9019]	568	LOGICAL :: lit
[4666]	569	REAL(wp) :: zlamb1, zlamb2, zlamb3
	570	REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
	571	REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
[6140]	572	REAL(wp) :: zeps = 1.e-20_wp
	573	REAL(wp) :: zerr
[9019]	574	REAL(wp), DIMENSION(jpi,jpj) :: zfrz
	575	REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas
	576	REAL(wp), DIMENSION(jpi,jpj) :: zfwflx, zhtflx, zhtflx_b
[4666]	577	!!---------------------------------------------------------------------
[9124]	578	!
[6140]	579	! coeficient for linearisation of potential tfreez
	580	! Crude approximation for pressure (but commonly used)
[9728]	581	IF ( l_useCT ) THEN ! linearisation from Jourdain et al. (2017)
	582	zlamb1 =-0.0564_wp
	583	zlamb2 = 0.0773_wp
	584	zlamb3 =-7.8633e-8 * grav * rau0
	585	ELSE ! linearisation from table 4 (Asay-Davis et al., 2015)
	586	zlamb1 =-0.0573_wp
	587	zlamb2 = 0.0832_wp
	588	zlamb3 =-7.53e-8 * grav * rau0
	589	ENDIF
[4666]	590	!
[6140]	591	! initialisation
	592	zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
	593	zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp
	594	zfwflx (:,:) = 0.0_wp
[4666]	595
[6140]	596	! compute ice shelf melting
	597	nit = 1 ; lit = .TRUE.
	598	DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine
	599	SELECT CASE ( nn_isfblk )
	600	CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
	601	! Calculate freezing temperature
	602	CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
[4666]	603
[6140]	604	! compute gammat every where (2d)
	605	CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
	606
	607	! compute upward heat flux zhtflx and upward water flux zwflx
	608	DO jj = 1, jpj
	609	DO ji = 1, jpi
	610	zhtflx(ji,jj) = zgammat(ji,jj)rcprau0*(ttbl(ji,jj)-zfrz(ji,jj))
[9935]	611	zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf
[6140]	612	END DO
	613	END DO
[4666]	614
[6140]	615	! Compute heat flux and upward fresh water flux
	616	qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
	617	fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
[5721]	618
[6140]	619	CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
	620	! compute gammat every where (2d)
	621	CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
[4666]	622
[6140]	623	! compute upward heat flux zhtflx and upward water flux zwflx
	624	! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
	625	DO jj = 1, jpj
	626	DO ji = 1, jpi
	627	! compute coeficient to solve the 2nd order equation
	628	zeps1 = rcprau0zgammat(ji,jj)
[9935]	629	zeps2 = rLfusisfrau0zgammas(ji,jj)
	630	zeps3 = rhoisfrcpisfrkappa/MAX(risfdep(ji,jj),zeps)
[6140]	631	zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
	632	zeps6 = zeps4-ttbl(ji,jj)
	633	zeps7 = zeps4-tsurf
	634	zaqe = zlamb1 * (zeps1 + zeps3)
	635	zaqer = 0.5_wp/MIN(zaqe,-zeps)
	636	zbqe = zeps1zeps6+zeps3zeps7-zeps2
	637	zcqe = zeps2*stbl(ji,jj)
	638	zdis = zbqezbqe-4.0_wpzaqe*zcqe
	639
	640	! Presumably zdis can never be negative because gammas is very small compared to gammat
	641	! compute s freeze
	642	zsfrz=(-zbqe-SQRT(zdis))*zaqer
	643	IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
	644
	645	! compute t freeze (eq. 22)
	646	zfrz(ji,jj)=zeps4+zlamb1*zsfrz
	647
	648	! zfwflx is upward water flux
	649	! zhtflx is upward heat flux (out of ocean)
	650	! compute the upward water and heat flux (eq. 28 and eq. 29)
	651	zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
	652	zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) )
	653	END DO
[4666]	654	END DO
	655
[6140]	656	! compute heat and water flux
	657	qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
	658	fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
[4666]	659
[6140]	660	END SELECT
[4666]	661
[6140]	662	! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
	663	IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE.
	664	ELSE
	665	! check total number of iteration
	666	IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
	667	ELSE ; nit = nit + 1
	668	END IF
[4666]	669
[6140]	670	! compute error between 2 iterations
	671	! if needed save gammat and compute zhtflx_b for next iteration
	672	zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
	673	IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
	674	ELSE ; zhtflx_b(:,:) = zhtflx(:,:)
	675	END IF
	676	END IF
	677	END DO
[5302]	678	!
[6140]	679	CALL iom_put('isfgammat', zgammat)
	680	CALL iom_put('isfgammas', zgammas)
	681	!
[4666]	682	END SUBROUTINE sbc_isf_cav
	683
[7646]	684
[6140]	685	SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
[4666]	686	!!----------------------------------------------------------------------
	687	!! ** Purpose : compute the coefficient echange for heat flux
	688	!!
	689	!! ** Method : gamma assume constant or depends of u* and stability
	690	!!
	691	!! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
	692	!! Jenkins et al., 2010, JPO, p2298-2312
	693	!!---------------------------------------------------------------------
[9124]	694	REAL(wp), DIMENSION(:,:), INTENT( out) :: pgt , pgs !
	695	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf !
[6140]	696	!
	697	INTEGER :: ji, jj ! loop index
[9124]	698	INTEGER :: ikt ! local integer
[4666]	699	REAL(wp) :: zdku, zdkv ! U, V shear
	700	REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number
	701	REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point
	702	REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness
	703	REAL(wp) :: zhmax ! limitation of mol
	704	REAL(wp) :: zetastar ! stability parameter
	705	REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence
	706	REAL(wp) :: zcoef ! temporary coef
[4946]	707	REAL(wp) :: zdep
[6140]	708	REAL(wp) :: zeps = 1.0e-20_wp
	709	REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant
	710	REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
[4946]	711	REAL(wp), DIMENSION(2) :: zts, zab
[9019]	712	REAL(wp), DIMENSION(jpi,jpj) :: zustar ! U, V at T point and friction velocity
[4666]	713	!!---------------------------------------------------------------------
	714	!
[6140]	715	SELECT CASE ( nn_gammablk )
	716	CASE ( 0 ) ! gamma is constant (specified in namelist)
	717	!! ISOMIP formulation (Hunter et al, 2006)
	718	pgt(:,:) = rn_gammat0
	719	pgs(:,:) = rn_gammas0
[4666]	720
[6140]	721	CASE ( 1 ) ! gamma is assume to be proportional to u*
	722	!! Jenkins et al., 2010, JPO, p2298-2312
	723	!! Adopted by Asay-Davis et al. (2015)
[9124]	724	!! compute ustar (eq. 24)
	725	!!gm NB use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
	726	!! zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
	727	!! or better : compute ustar in zdfdrg and use it here as well as in TKE, GLS and Co
	728	!!
	729	!! ===>>>> GM to be done this chrismas
	730	!!
	731	!!gm end
[9019]	732	zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
[4666]	733
[6140]	734	!! Compute gammats
	735	pgt(:,:) = zustar(:,:) * rn_gammat0
	736	pgs(:,:) = zustar(:,:) * rn_gammas0
	737
	738	CASE ( 2 ) ! gamma depends of stability of boundary layer
	739	!! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
	740	!! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
	741	!! compute ustar
[9019]	742	zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
[4666]	743
[6140]	744	!! compute Pr and Sc number (can be improved)
	745	zPr = 13.8_wp
	746	zSc = 2432.0_wp
[4666]	747
[6140]	748	!! compute gamma mole
	749	zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
	750	zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
[4666]	751
[6140]	752	!! compute gamma
[9019]	753	DO ji = 2, jpi
	754	DO jj = 2, jpj
[6140]	755	ikt = mikt(ji,jj)
[4666]	756
[9019]	757	IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think
[6140]	758	pgt = rn_gammat0
	759	pgs = rn_gammas0
	760	ELSE
	761	!! compute Rc number (as done in zdfric.F90)
[9019]	762	!!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation
	763	!!gm moreover, use Max(rn2,0) to take care of static instabilities....
[10536]	764	zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
[6140]	765	! ! shear of horizontal velocity
	766	zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) &
	767	& -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) )
	768	zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) &
	769	& -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) )
	770	! ! richardson number (minimum value set to zero)
	771	zRc = rn2(ji,jj,ikt+1) / MAX( zdkuzdku + zdkvzdkv, zeps )
	772
	773	!! compute bouyancy
	774	zts(jp_tem) = ttbl(ji,jj)
	775	zts(jp_sal) = stbl(ji,jj)
	776	zdep = gdepw_n(ji,jj,ikt)
[4946]	777	!
[6140]	778	CALL eos_rab( zts, zdep, zab )
	779	!
	780	!! compute length scale
	781	zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!
[4666]	782
[6140]	783	!! compute Monin Obukov Length
	784	! Maximum boundary layer depth
	785	zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
	786	! Compute Monin obukhov length scale at the surface and Ekman depth:
	787	zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
	788	zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
[4666]	789
[6140]	790	!! compute eta* (stability parameter)
[7646]	791	zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
[4666]	792
[6140]	793	!! compute the sublayer thickness
	794	zhnu = 5 * znu / zustar(ji,jj)
[4666]	795
[6140]	796	!! compute gamma turb
[7646]	797	zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
[6140]	798	& + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
	799
	800	!! compute gammats
	801	pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
	802	pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
[4666]	803	END IF
[6140]	804	END DO
	805	END DO
[10425]	806	CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.)
[6140]	807	END SELECT
[5836]	808	!
[6140]	809	END SUBROUTINE sbc_isf_gammats
[4666]	810
[9019]	811
[6140]	812	SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
[4666]	813	!!----------------------------------------------------------------------
	814	!! * SUBROUTINE sbc_isf_tbl *
	815	!!
[6140]	816	!! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
[4666]	817	!!
	818	!!----------------------------------------------------------------------
[9019]	819	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvarin
	820	REAL(wp), DIMENSION(:,:) , INTENT( out) :: pvarout
	821	CHARACTER(len=1), INTENT(in ) :: cd_ptin ! point of variable in/out
[6140]	822	!
[9019]	823	INTEGER :: ji, jj, jk ! loop index
	824	INTEGER :: ikt, ikb ! top and bottom index of the tbl
	825	REAL(wp) :: ze3, zhk
	826	REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
[6140]	827	!!----------------------------------------------------------------------
[4666]	828
[6140]	829	! initialisation
	830	pvarout(:,:)=0._wp
	831
	832	SELECT CASE ( cd_ptin )
	833	CASE ( 'U' ) ! compute U in the top boundary layer at T- point
	834	DO jj = 1,jpj
	835	DO ji = 1,jpi
	836	ikt = miku(ji,jj) ; ikb = miku(ji,jj)
	837	! thickness of boundary layer at least the top level thickness
[9019]	838	zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
[4666]	839
[6140]	840	! determine the deepest level influenced by the boundary layer
	841	DO jk = ikt+1, mbku(ji,jj)
	842	IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
	843	END DO
	844	zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
[4666]	845
[6140]	846	! level fully include in the ice shelf boundary layer
	847	DO jk = ikt, ikb - 1
	848	ze3 = e3u_n(ji,jj,jk)
	849	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
	850	END DO
[4666]	851
[6140]	852	! level partially include in ice shelf boundary layer
	853	zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
	854	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
	855	END DO
	856	END DO
[9019]	857	DO jj = 2, jpj
	858	DO ji = 2, jpi
	859	!!gm a wet-point only average should be used here !!!
[6140]	860	pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
	861	END DO
	862	END DO
[10425]	863	CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
[6140]	864
	865	CASE ( 'V' ) ! compute V in the top boundary layer at T- point
	866	DO jj = 1,jpj
	867	DO ji = 1,jpi
	868	ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
	869	! thickness of boundary layer at least the top level thickness
	870	zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
[4666]	871
[6140]	872	! determine the deepest level influenced by the boundary layer
	873	DO jk = ikt+1, mbkv(ji,jj)
	874	IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
	875	END DO
	876	zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
[4666]	877
[6140]	878	! level fully include in the ice shelf boundary layer
	879	DO jk = ikt, ikb - 1
	880	ze3 = e3v_n(ji,jj,jk)
	881	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
	882	END DO
[4666]	883
[6140]	884	! level partially include in ice shelf boundary layer
	885	zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
	886	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
	887	END DO
	888	END DO
[9019]	889	DO jj = 2, jpj
	890	DO ji = 2, jpi
	891	!!gm a wet-point only average should be used here !!!
[6140]	892	pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
	893	END DO
	894	END DO
[10425]	895	CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
[6140]	896
	897	CASE ( 'T' ) ! compute T in the top boundary layer at T- point
	898	DO jj = 1,jpj
	899	DO ji = 1,jpi
	900	ikt = misfkt(ji,jj)
	901	ikb = misfkb(ji,jj)
	902
[4666]	903	! level fully include in the ice shelf boundary layer
	904	DO jk = ikt, ikb - 1
[6140]	905	ze3 = e3t_n(ji,jj,jk)
	906	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
[4666]	907	END DO
	908
	909	! level partially include in ice shelf boundary layer
[6140]	910	zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
	911	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
	912	END DO
[4666]	913	END DO
[6140]	914	END SELECT
[9019]	915	!
[6140]	916	! mask mean tbl value
	917	pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
[5836]	918	!
[4666]	919	END SUBROUTINE sbc_isf_tbl
	920
	921
	922	SUBROUTINE sbc_isf_div( phdivn )
	923	!!----------------------------------------------------------------------
	924	!! * SUBROUTINE sbc_isf_div *
	925	!!
	926	!! ** Purpose : update the horizontal divergence with the runoff inflow
	927	!!
	928	!! ** Method :
	929	!! CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
	930	!! divergence and expressed in m/s
	931	!!
	932	!! ** Action : phdivn decreased by the runoff inflow
	933	!!----------------------------------------------------------------------
[6140]	934	REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence
	935	!
[4946]	936	INTEGER :: ji, jj, jk ! dummy loop indices
	937	INTEGER :: ikt, ikb
[6140]	938	REAL(wp) :: zhk
	939	REAL(wp) :: zfact ! local scalar
[4666]	940	!!----------------------------------------------------------------------
	941	!
	942	zfact = 0.5_wp
	943	!
[6140]	944	IF(.NOT.ln_linssh ) THEN ! need to re compute level distribution of isf fresh water
[4666]	945	DO jj = 1,jpj
	946	DO ji = 1,jpi
	947	ikt = misfkt(ji,jj)
	948	ikb = misfkt(ji,jj)
	949	! thickness of boundary layer at least the top level thickness
[6140]	950	rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
[4666]	951
	952	! determine the deepest level influenced by the boundary layer
	953	DO jk = ikt, mbkt(ji,jj)
[6140]	954	IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
[4666]	955	END DO
[6140]	956	rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
[4666]	957	misfkb(ji,jj) = ikb ! last wet level of the tbl
	958	r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
[4726]	959
[6140]	960	zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
	961	ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer
[4666]	962	END DO
	963	END DO
[6140]	964	END IF
[4666]	965	!
[6140]	966	!== ice shelf melting distributed over several levels ==!
[4666]	967	DO jj = 1,jpj
	968	DO ji = 1,jpi
	969	ikt = misfkt(ji,jj)
	970	ikb = misfkb(ji,jj)
	971	! level fully include in the ice shelf boundary layer
	972	DO jk = ikt, ikb - 1
[4946]	973	phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
[6140]	974	& * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
[4666]	975	END DO
	976	! level partially include in ice shelf boundary layer
[4946]	977	phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
[6140]	978	& + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)
[4666]	979	END DO
	980	END DO
	981	!
	982	END SUBROUTINE sbc_isf_div
[9019]	983
[6140]	984	!!======================================================================
[4666]	985	END MODULE sbcisf

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

Download in other formats: