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.

sbcblk.F90 in NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90 @ 12015

Last change on this file since 12015 was 12015, checked in by gsamson, 4 years ago
dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK
Property svn:keywords set to `Id`
File size: 76.6 KB

Rev	Line
[6723]	1	MODULE sbcblk
	2	!!======================================================================
	3	!! * MODULE sbcblk *
	4	!! Ocean forcing: momentum, heat and freshwater flux formulation
	5	!! Aerodynamic Bulk Formulas
	6	!! SUCCESSOR OF "sbcblk_core"
	7	!!=====================================================================
[7163]	8	!! History : 1.0 ! 2004-08 (U. Schweckendiek) Original CORE code
	9	!! 2.0 ! 2005-04 (L. Brodeau, A.M. Treguier) improved CORE bulk and its user interface
	10	!! 3.0 ! 2006-06 (G. Madec) sbc rewritting
	11	!! - ! 2006-12 (L. Brodeau) Original code for turb_core
[6723]	12	!! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put
	13	!! 3.3 ! 2010-10 (S. Masson) add diurnal cycle
[7163]	14	!! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE
	15	!! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk
	16	!! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
[12015]	17	!! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
[7163]	18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
[12015]	19	!! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle)
[11275]	20	!! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE)
[6723]	21	!!----------------------------------------------------------------------
	22
	23	!!----------------------------------------------------------------------
[7163]	24	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
	25	!! sbc_blk : bulk formulation as ocean surface boundary condition
[12015]	26	!! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE)
	27	!! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE)
	28	!! sea-ice case only :
[11304]	29	!! blk_ice_1 : provide the air-ice stress
	30	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
[10534]	31	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	32	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
[12015]	33	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[6723]	34	!!----------------------------------------------------------------------
	35	USE oce ! ocean dynamics and tracers
	36	USE dom_oce ! ocean space and time domain
	37	USE phycst ! physical constants
	38	USE fldread ! read input fields
	39	USE sbc_oce ! Surface boundary condition: ocean fields
	40	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
	41	USE sbcdcy ! surface boundary condition: diurnal cycle
	42	USE sbcwave , ONLY : cdn_wave ! wave module
	43	USE sbc_ice ! Surface boundary condition: ice fields
	44	USE lib_fortran ! to use key_nosignedzero
[9570]	45	#if defined key_si3
[11360]	46	USE ice , ONLY : jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif
[9019]	47	USE icethd_dh ! for CALL ice_thd_snwblow
[6723]	48	#endif
[12015]	49	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)
	50	USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
	51	USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
	52	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1)
[6727]	53	!
[6723]	54	USE iom ! I/O manager library
	55	USE in_out_manager ! I/O manager
	56	USE lib_mpp ! distribued memory computing library
	57	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	58	USE prtctl ! Print control
	59
[12015]	60	USE sbcblk_phy !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc...
	61
	62
[6723]	63	IMPLICIT NONE
	64	PRIVATE
	65
[7163]	66	PUBLIC sbc_blk_init ! called in sbcmod
	67	PUBLIC sbc_blk ! called in sbcmod
[11275]	68	PUBLIC blk_oce_1 ! called in sbcabl
	69	PUBLIC blk_oce_2 ! called in sbcabl
[9570]	70	#if defined key_si3
[11304]	71	PUBLIC blk_ice_1 ! routine called in icesbc
	72	PUBLIC blk_ice_2 ! routine called in icesbc
[10535]	73	PUBLIC blk_ice_qcn ! routine called in icesbc
[12015]	74	#endif
[6723]	75
[12015]	76	INTEGER , PUBLIC :: jpfld ! maximum number of files to read
	77	INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
	78	INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
	79	INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
	80	INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
	81	INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
	82	INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
	83	INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
	84	INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
	85	INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa)
	86	INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
	87	INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
[11275]	88
[12015]	89	TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read)
[6723]	90
[12015]	91	! !!* Namelist namsbc_blk : bulk parameters
	92	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
	93	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
	94	LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013)
	95	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1)
[11309]	96	!
[12015]	97	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
	98	REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation
	99	REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress
	100	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
	101	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
	102	!
	103	LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012)
	104	LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
	105	!
[11360]	106	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice
	107	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme)
	108	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
[11309]	109
[12015]	110	LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
	111	LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
	112	LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
	113	LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
	114	LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB
	115	!
	116	INTEGER :: nhumi ! choice of the bulk algorithm
	117	! ! associated indices:
	118	INTEGER, PARAMETER :: np_humi_sph = 1
	119	INTEGER, PARAMETER :: np_humi_dpt = 2
	120	INTEGER, PARAMETER :: np_humi_rlh = 3
[6723]	121
[12015]	122	INTEGER :: nblk ! choice of the bulk algorithm
	123	! ! associated indices:
	124	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
	125	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
	126	INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013)
	127	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1)
[6723]	128
	129	!! * Substitutions
	130	# include "vectopt_loop_substitute.h90"
	131	!!----------------------------------------------------------------------
[9598]	132	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	133	!! $Id$
[10068]	134	!! Software governed by the CeCILL license (see ./LICENSE)
[6723]	135	!!----------------------------------------------------------------------
	136	CONTAINS
	137
[7355]	138	INTEGER FUNCTION sbc_blk_alloc()
	139	!!-------------------------------------------------------------------
	140	!! * ROUTINE sbc_blk_alloc *
	141	!!-------------------------------------------------------------------
[11360]	142	ALLOCATE( t_zu(jpi,jpj), q_zu(jpi,jpj), &
	143	& Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), &
	144	& Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc )
[7355]	145	!
[10425]	146	CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
	147	IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
[7355]	148	END FUNCTION sbc_blk_alloc
	149
[9019]	150
[7163]	151	SUBROUTINE sbc_blk_init
	152	!!---------------------------------------------------------------------
	153	!! * ROUTINE sbc_blk_init *
	154	!!
	155	!! ** Purpose : choose and initialize a bulk formulae formulation
	156	!!
[12015]	157	!! ** Method :
[7163]	158	!!
	159	!!----------------------------------------------------------------------
[11309]	160	INTEGER :: jfpr ! dummy loop indice and argument
[7163]	161	INTEGER :: ios, ierror, ioptio ! Local integer
	162	!!
	163	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
[12015]	164	TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read
[7163]	165	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
	166	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
[11275]	167	TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " "
[7163]	168	NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields
[11275]	169	& sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, &
[12015]	170	& ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm
	171	& cn_dir , rn_zqt, rn_zu, &
	172	& rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, &
	173	& ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB
[7163]	174	!!---------------------------------------------------------------------
	175	!
[7355]	176	! ! allocate sbc_blk_core array
	177	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
	178	!
[12015]	179	! !** read bulk namelist
[7163]	180	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
	181	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
[11348]	182	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
[7163]	183	!
	184	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
	185	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
[11348]	186	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
[7163]	187	!
	188	IF(lwm) WRITE( numond, namsbc_blk )
	189	!
	190	! !** initialization of the chosen bulk formulae (+ check)
	191	! !* select the bulk chosen in the namelist and check the choice
[12015]	192	ioptio = 0
	193	IF( ln_NCAR ) THEN
	194	nblk = np_NCAR ; ioptio = ioptio + 1
	195	ENDIF
	196	IF( ln_COARE_3p0 ) THEN
	197	nblk = np_COARE_3p0 ; ioptio = ioptio + 1
	198	ENDIF
	199	IF( ln_COARE_3p6 ) THEN
	200	nblk = np_COARE_3p6 ; ioptio = ioptio + 1
	201	ENDIF
	202	IF( ln_ECMWF ) THEN
	203	nblk = np_ECMWF ; ioptio = ioptio + 1
	204	ENDIF
[7163]	205	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
[12015]	206
	207	! !** initialization of the cool-skin / warm-layer parametrization
	208	IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) &
	209	& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
[7163]	210	!
[12015]	211	ioptio = 0
	212	IF( ln_humi_sph ) THEN
	213	nhumi = np_humi_sph ; ioptio = ioptio + 1
	214	ENDIF
	215	IF( ln_humi_dpt ) THEN
	216	nhumi = np_humi_dpt ; ioptio = ioptio + 1
	217	ENDIF
	218	IF( ln_humi_rlh ) THEN
	219	nhumi = np_humi_rlh ; ioptio = ioptio + 1
	220	ENDIF
	221	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
	222	!
[7163]	223	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
[11267]	224	IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
[12015]	225	IF( sn_qsr%ln_tint ) THEN
[7163]	226	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
	227	& ' ==> We force time interpolation = .false. for qsr' )
	228	sn_qsr%ln_tint = .false.
	229	ENDIF
	230	ENDIF
	231	! !* set the bulk structure
	232	! !- store namelist information in an array
[11586]	233	IF( ln_blk ) jpfld = 9
	234	IF( ln_abl ) jpfld = 11
	235	ALLOCATE( slf_i(jpfld) )
	236	!
[7163]	237	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
	238	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
	239	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
	240	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
[11304]	241	slf_i(jp_slp ) = sn_slp
[11586]	242	IF( ln_abl ) THEN
	243	slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj
	244	END IF
[7163]	245	!
	246	! !- allocate the bulk structure
[11275]	247	ALLOCATE( sf(jpfld), STAT=ierror )
[7163]	248	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
	249	!
[11309]	250	DO jfpr= 1, jpfld
[11275]	251	!
[11309]	252	IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero)
	253	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
	254	sf(jfpr)%fnow(:,:,1) = 0._wp
	255	ELSE !-- used field --!
	256	IF( ln_abl .AND. &
	257	& ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. &
	258	& jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN ! ABL: some fields are 3D input
	259	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) )
	260	IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) )
	261	ELSE ! others or Bulk fields are 2D fiels
	262	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
	263	IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) )
[11275]	264	ENDIF
	265	!
[11309]	266	IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) &
[11275]	267	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
	268	& ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
	269	ENDIF
	270	END DO
[12015]	271	! !- fill the bulk structure with namelist informations
	272	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
[11275]	273	!
[12015]	274	IF( ln_wave ) THEN
	275	!Activated wave module but neither drag nor stokes drift activated
	276	IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN
	277	CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
	278	!drag coefficient read from wave model definable only with mfs bulk formulae and core
	279	ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN
	280	CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
	281	ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN
	282	CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
	283	ENDIF
[11309]	284	ELSE
[12015]	285	IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &
	286	& CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &
	287	& 'with drag coefficient (ln_cdgw =T) ' , &
	288	& 'or Stokes Drift (ln_sdw=T) ' , &
	289	& 'or ocean stress modification due to waves (ln_tauwoc=T) ', &
	290	& 'or Stokes-Coriolis term (ln_stcori=T)' )
	291	ENDIF
[7431]	292	!
[11309]	293	IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient
[11275]	294	rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level
	295	rn_zu = ght_abl(2)
	296	IF(lwp) WRITE(numout,*)
[11360]	297	IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
[11275]	298	ENDIF
[11360]	299	!
	300	! set transfer coefficients to default sea-ice values
[12015]	301	Cd_ice(:,:) = rCd_ice
	302	Ch_ice(:,:) = rCd_ice
	303	Ce_ice(:,:) = rCd_ice
[11360]	304	!
[7163]	305	IF(lwp) THEN !** Control print
	306	!
[12015]	307	WRITE(numout,) ! namelist
[7163]	308	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
	309	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
	310	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
[12015]	311	WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
	312	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF
[7163]	313	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
	314	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
	315	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
	316	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
	317	WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
	318	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
[9019]	319	WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12
	320	WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15
[7163]	321	!
	322	WRITE(numout,*)
	323	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
[9190]	324	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
	325	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
[12015]	326	CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
	327	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)'
[7163]	328	END SELECT
	329	!
[12015]	330	WRITE(numout,*)
	331	WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs
	332	WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl
	333	!
	334	WRITE(numout,*)
	335	SELECT CASE( nhumi ) !* Print the choice of air humidity
	336	CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]'
	337	CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]'
	338	CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]'
	339	END SELECT
	340	!
[7163]	341	ENDIF
	342	!
	343	END SUBROUTINE sbc_blk_init
	344
	345
[6723]	346	SUBROUTINE sbc_blk( kt )
	347	!!---------------------------------------------------------------------
	348	!! * ROUTINE sbc_blk *
	349	!!
	350	!! ** Purpose : provide at each time step the surface ocean fluxes
[9019]	351	!! (momentum, heat, freshwater and runoff)
[6723]	352	!!
	353	!! ** Method : (1) READ each fluxes in NetCDF files:
	354	!! the 10m wind velocity (i-component) (m/s) at T-point
	355	!! the 10m wind velocity (j-component) (m/s) at T-point
[11275]	356	!! the 10m or 2m specific humidity (kg/kg)
[6723]	357	!! the solar heat (W/m2)
	358	!! the Long wave (W/m2)
	359	!! the 10m or 2m air temperature (Kelvin)
	360	!! the total precipitation (rain+snow) (Kg/m2/s)
[11275]	361	!! the snow (solid precipitation) (kg/m2/s)
	362	!! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)
	363	!! (2) CALL blk_oce_1 and blk_oce_2
[6723]	364	!!
	365	!! C A U T I O N : never mask the surface stress fields
	366	!! the stress is assumed to be in the (i,j) mesh referential
	367	!!
	368	!! ** Action : defined at each time-step at the air-sea interface
	369	!! - utau, vtau i- and j-component of the wind stress
	370	!! - taum wind stress module at T-point
	371	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
	372	!! - qns, qsr non-solar and solar heat fluxes
	373	!! - emp upward mass flux (evapo. - precip.)
	374	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
	375	!!
	376	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
	377	!! Brodeau et al. Ocean Modelling 2010
	378	!!----------------------------------------------------------------------
	379	INTEGER, INTENT(in) :: kt ! ocean time step
[11275]	380	!!
	381	REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp
[6723]	382	!!---------------------------------------------------------------------
	383	!
	384	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
[7163]	385	!
[12015]	386	IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1) ! no previous estimate of skin temperature => using bulk SST
	387	!
	388	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[11275]	389	CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in
	390	& sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in
	391	& sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in
[12015]	392	& sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs)
[11275]	393	& zssq, zcd_du, zsen, zevp ) ! =>> out
	394
[11360]	395	CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in
[11275]	396	& sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in
	397	& sf(jp_snow)%fnow(:,:,1), sst_m, & ! <<= in
	398	& zsen, zevp ) ! <=> in out
	399	ENDIF
[12015]	400	!
[6723]	401	#if defined key_cice
[12015]	402	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[7753]	403	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
[12015]	404	IF( ln_dm2dc ) THEN
	405	qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
	406	ELSE
	407	qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
	408	ENDIF
[7753]	409	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
[12015]	410
	411	SELECT CASE( nhumi )
	412	CASE( np_humi_sph )
	413	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
	414	CASE( np_humi_dpt )
	415	qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
	416	CASE( np_humi_rlh )
	417	qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file
	418	END SELECT
	419
[7753]	420	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
	421	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
	422	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
	423	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
[6723]	424	ENDIF
	425	#endif
	426	!
	427	END SUBROUTINE sbc_blk
	428
	429
[12015]	430	SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp
	431	& pslp , pst , pu , pv, & ! inp
	432	& pqsr , pqlw , & ! inp
	433	& pssq , pcd_du, psen , pevp ) ! out
[6723]	434	!!---------------------------------------------------------------------
[11275]	435	!! * ROUTINE blk_oce_1 *
[6723]	436	!!
[12015]	437	!! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes
	438	!! if ln_abl=T, computes Cd x \|U\|, Ch x \|U\|, Ce x \|U\| for ABL integration
[6723]	439	!!
[12015]	440	!! ** Method : bulk formulae using atmospheric fields from :
	441	!! if ln_blk=T, atmospheric fields read in sbc_read
	442	!! if ln_abl=T, the ABL model at previous time-step
[6723]	443	!!
[11275]	444	!! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg)
	445	!! - pcd_du : Cd x \|dU\| at T-points (m/s)
	446	!! - psen : Ch x \|dU\| at T-points (m/s)
	447	!! - pevp : Ce x \|dU\| at T-points (m/s)
[6723]	448	!!---------------------------------------------------------------------
[11275]	449	INTEGER , INTENT(in ) :: kt ! time step index
	450	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s]
	451	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s]
	452	REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg]
	453	REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin]
	454	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa]
	455	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celcius]
	456	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
	457	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
[12015]	458	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr !
	459	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw !
[11275]	460	REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg]
	461	REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x \|dU\| at T-points [m/s]
	462	REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x \|dU\| at T-points [m/s]
	463	REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x \|dU\| at T-points [m/s]
[6723]	464	!
	465	INTEGER :: ji, jj ! dummy loop indices
	466	REAL(wp) :: zztmp ! local variable
[9019]	467	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
	468	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
	469	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
	470	REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]
[12015]	471	REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg]
[11360]	472	REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
	473	REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
	474	REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean
[12015]	475	REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux
	476	REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2
[6723]	477	!!---------------------------------------------------------------------
	478	!
[7753]	479	! local scalars ( place there for vector optimisation purposes)
	480	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
[6723]	481
	482	! ----------------------------------------------------------------------------- !
	483	! 0 Wind components and module at T-point relative to the moving ocean !
	484	! ----------------------------------------------------------------------------- !
	485
[7753]	486	! ... components ( U10m - U_oce ) at T-point (unmasked)
[11275]	487	#if defined key_cyclone
[7753]	488	zwnd_i(:,:) = 0._wp
	489	zwnd_j(:,:) = 0._wp
[6723]	490	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
	491	DO jj = 2, jpjm1
	492	DO ji = fs_2, fs_jpim1 ! vect. opt.
[11275]	493	pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
	494	pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
[6723]	495	END DO
	496	END DO
	497	#endif
	498	DO jj = 2, jpjm1
	499	DO ji = fs_2, fs_jpim1 ! vect. opt.
[11275]	500	zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
	501	zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
[6723]	502	END DO
	503	END DO
[10425]	504	CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
[6723]	505	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
[7753]	506	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
	507	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
[6723]	508
	509	! ----------------------------------------------------------------------------- !
[12015]	510	! I Solar FLUX !
[6723]	511	! ----------------------------------------------------------------------------- !
	512
[12015]	513	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
	514	zztmp = 1. - albo
	515	IF( ln_dm2dc ) THEN
	516	qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
	517	ELSE
	518	qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
	519	ENDIF
	520
	521
	522	! ----------------------------------------------------------------------------- !
	523	! II Turbulent FLUXES !
	524	! ----------------------------------------------------------------------------- !
	525
	526	! specific humidity at SST
	527	pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) )
	528
	529	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
	530	zztmp1(:,:) = zst(:,:)
	531	zztmp2(:,:) = pssq(:,:)
	532	ENDIF
	533
	534	! specific humidity of air at "rn_zqt" m above the sea
	535	SELECT CASE( nhumi )
	536	CASE( np_humi_sph )
	537	zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity!
	538	CASE( np_humi_dpt )
	539	!IF(lwp) WRITE(numout,) ' ** blk_oce => computing q_air out of d_air and slp !' !LBrm
	540	zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) )
	541	CASE( np_humi_rlh )
	542	!IF(lwp) WRITE(numout,) ' ** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
	543	zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file
	544	END SELECT
[11275]	545	!
[12015]	546	! potential temperature of air at "rn_zqt" m above the sea
[11275]	547	IF( ln_abl ) THEN
	548	ztpot = ptair(:,:)
	549	ELSE
	550	! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
	551	! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
	552	! (since reanalysis products provide T at z, not theta !)
[12015]	553	ztpot = ptair(:,:) + gamma_moist( ptair(:,:), zqair(:,:) ) * rn_zqt
[6723]	554	ENDIF
	555
	556
	557
[12015]	558	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
[7355]	559
[12015]	560	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
	561
	562	CASE( np_COARE_3p0 )
	563	CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	564	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
	565	& Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
	566
	567	CASE( np_COARE_3p6 )
	568	CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	569	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
	570	& Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
	571
	572	CASE( np_ECMWF )
	573	CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	574	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
	575	& Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
	576
	577	CASE DEFAULT
	578	CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' )
	579	END SELECT
	580
	581	!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq:
	582	WHERE ( fr_i < 0.001_wp )
	583	! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!!
	584	zst(:,:) = zst(:,:)*tmask(:,:,1)
	585	pssq(:,:) = pssq(:,:)*tmask(:,:,1)
	586	ELSEWHERE
	587	! we forget about the update...
	588	zst(:,:) = zztmp1(:,:) !LB: using what we backed up before skin-algo
	589	pssq(:,:) = zztmp2(:,:) !LB: " " "
	590	END WHERE
	591
	592	!LB: Update of tsk, the "official" array for skin temperature
	593	tsk(:,:) = zst(:,:)
	594
	595	ELSE !IF( ln_skin_cs .OR. ln_skin_wl )
	596
	597	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
	598	!
	599	CASE( np_NCAR )
	600	CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, &
	601	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
	602
	603	CASE( np_COARE_3p0 )
	604	CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	605	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
	606
	607	CASE( np_COARE_3p6 )
	608	CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	609	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
	610
	611	CASE( np_ECMWF )
	612	CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
	613	& zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
	614
	615	CASE DEFAULT
	616	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
	617	END SELECT
	618
	619	ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl )
	620
	621	!! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef.
	622	!! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef.
	623
	624	IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
	625	!! If zu == zt, then ensuring once for all that:
	626	t_zu(:,:) = ztpot(:,:)
	627	q_zu(:,:) = zqair(:,:)
	628	ENDIF
	629
	630
	631	! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90
	632	! -------------------------------------------------------------
	633
	634	IF (ln_abl) THEN
	635
	636	CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
	637	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
	638	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
	639	& pcd_du(:,:), psen(:,:), pevp(:,:), &
	640	& prhoa=rhoa(:,:) )
	641
	642	ELSE
	643
	644	CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
	645	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
	646	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
	647	& taum(:,:), psen(:,:), zqla(:,:), &
	648	& pEvap=pevp(:,:), prhoa=rhoa(:,:) )
	649
	650	zqla(:,:) = zqla(:,:) * tmask(:,:,1)
	651	psen(:,:) = psen(:,:) * tmask(:,:,1)
	652	taum(:,:) = taum(:,:) * tmask(:,:,1)
	653	pevp(:,:) = pevp(:,:) * tmask(:,:,1)
[11275]	654
[12015]	655	! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array...
	656	zcd_oce = 0._wp
	657	WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm
	658	zwnd_i = zcd_oce * zwnd_i
	659	zwnd_j = zcd_oce * zwnd_j
[6723]	660
[11275]	661	CALL iom_put( "taum_oce", taum ) ! output wind stress module
[12015]	662
[11275]	663	! ... utau, vtau at U- and V_points, resp.
	664	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
	665	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
	666	DO jj = 1, jpjm1
	667	DO ji = 1, fs_jpim1
	668	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
	669	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
	670	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
	671	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
	672	END DO
	673	END DO
	674	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
	675	IF(ln_ctl) THEN
	676	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
	677	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
	678	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
	679	ENDIF
	680	!
	681	ENDIF
	682	!
	683	IF(ln_ctl) THEN
	684	CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' )
	685	CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' )
	686	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' )
	687	ENDIF
	688	!
[12015]	689	END SUBROUTINE blk_oce_1
[6723]	690
	691
[11360]	692	SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in
	693	& psnow, pst , psen, pevp ) ! <<= in
[11275]	694	!!---------------------------------------------------------------------
	695	!! * ROUTINE blk_oce_2 *
	696	!!
	697	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
	698	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
	699	!! atmospheric variables (from ABL or external data)
	700	!!
	701	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
	702	!! - vtau : j-component of the stress at V-point (N/m2)
	703	!! - taum : Wind stress module at T-point (N/m2)
	704	!! - wndm : Wind speed module at T-point (m/s)
	705	!! - qsr : Solar heat flux over the ocean (W/m2)
	706	!! - qns : Non Solar heat flux over the ocean (W/m2)
	707	!! - emp : evaporation minus precipitation (kg/m2/s)
	708	!!---------------------------------------------------------------------
	709	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair
	710	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr
	711	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw
	712	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
	713	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
	714	REAL(wp), INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius]
	715	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
	716	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
	717	!
	718	INTEGER :: ji, jj ! dummy loop indices
	719	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
	720	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes
	721	REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation
	722	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
	723	!!---------------------------------------------------------------------
	724	!
	725	! local scalars ( place there for vector optimisation purposes)
	726	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
[6723]	727
[12015]	728
[11275]	729	! ----------------------------------------------------------------------------- !
[12015]	730	! III Net longwave radiative FLUX !
[11275]	731	! ----------------------------------------------------------------------------- !
[6723]	732
[12015]	733	!! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
	734	!! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
	735	zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefanzst(:,:)zst(:,:)zst(:,:)zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux
[6723]	736
[11275]	737	! Turbulent fluxes over ocean
	738	! -----------------------------
[6723]	739
[12015]	740	zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
[11275]	741
[6723]	742	IF(ln_ctl) THEN
[11275]	743	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' )
	744	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
[12015]	745
[6723]	746	ENDIF
	747
	748	! ----------------------------------------------------------------------------- !
[12015]	749	! IV Total FLUXES !
[6723]	750	! ----------------------------------------------------------------------------- !
	751	!
[12015]	752	emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.)
	753	& - pprec(:,:) * rn_pfac ) * tmask(:,:,1)
[7753]	754	!
[12015]	755	qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar
	756	& - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
	757	& - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !?
	758	& + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair
	759	& * ( ptair(:,:) - rt0 ) * rcp &
	760	& + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
	761	& * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
	762	qns(:,:) = qns(:,:) * tmask(:,:,1)
[7753]	763	!
[9570]	764	#if defined key_si3
[12015]	765	qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3)
[7753]	766	qsr_oce(:,:) = qsr(:,:)
[6723]	767	#endif
	768	!
[12015]	769	CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB
	770	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
	771	CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean
	772	CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean
	773	CALL iom_put( "evap_oce" , pevp ) ! evaporation
	774	CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean
	775	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
	776	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
	777	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
	778	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
	779	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
	780	CALL iom_put( 'snowpre', sprecip ) ! Snow
	781	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
	782	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
	783	CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius
	784	CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference...
[6723]	785	ENDIF
	786	!
	787	IF(ln_ctl) THEN
[11275]	788	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
	789	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
	790	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
[11360]	791	ENDIF
	792	!
[11275]	793	END SUBROUTINE blk_oce_2
[6723]	794
[11275]	795
[9570]	796	#if defined key_si3
[9019]	797	!!----------------------------------------------------------------------
[9570]	798	!! 'key_si3' SI3 sea-ice model
[9019]	799	!!----------------------------------------------------------------------
[11304]	800	!! blk_ice_1 : provide the air-ice stress
	801	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
[10534]	802	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	803	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
[12015]	804	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[9019]	805	!!----------------------------------------------------------------------
	806
[11360]	807	SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs
	808	& putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs
[6723]	809	!!---------------------------------------------------------------------
[11304]	810	!! * ROUTINE blk_ice_1 *
[6723]	811	!!
	812	!! ** Purpose : provide the surface boundary condition over sea-ice
	813	!!
	814	!! ** Method : compute momentum using bulk formulation
	815	!! formulea, ice variables and read atmospheric fields.
	816	!! NB: ice drag coefficient is assumed to be a constant
	817	!!---------------------------------------------------------------------
[11304]	818	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
	819	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
	820	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
	821	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
	822	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s]
	823	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
	824	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
[11360]	825	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
[11304]	826	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk
	827	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk
	828	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl
	829	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl
	830	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl
	831	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
	832	!
[6723]	833	INTEGER :: ji, jj ! dummy loop indices
[11360]	834	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
	835	REAL(wp) :: zootm_su ! sea-ice surface mean temperature
	836	REAL(wp) :: zztmp1, zztmp2 ! temporary arrays
	837	REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau)
[6723]	838	!!---------------------------------------------------------------------
[12015]	839	!
[6723]	840
[9019]	841	! ------------------------------------------------------------ !
	842	! Wind module relative to the moving ice ( U10m - U_ice ) !
	843	! ------------------------------------------------------------ !
[9767]	844	! C-grid ice dynamics : U & V-points (same as ocean)
	845	DO jj = 2, jpjm1
	846	DO ji = fs_2, fs_jpim1 ! vect. opt.
[11304]	847	zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) )
	848	zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
[9767]	849	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
[9019]	850	END DO
[9767]	851	END DO
[10425]	852	CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
[9767]	853	!
[9019]	854	! Make ice-atm. drag dependent on ice concentration
	855	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
[11360]	856	CALL Cdn10_Lupkes2012( Cd_ice )
	857	Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical
	858	Ce_ice(:,:) = Cd_ice(:,:)
[9019]	859	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
[12015]	860	CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice )
[11360]	861	Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical
[7355]	862	ENDIF
	863
[12015]	864	!! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef.
	865	!! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef.
[9019]	866
[6723]	867	! local scalars ( place there for vector optimisation purposes)
[12015]	868	IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)
[11360]	869	zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
[11304]	870
[11334]	871	IF( ln_blk ) THEN
[11304]	872	! ------------------------------------------------------------ !
	873	! Wind stress relative to the moving ice ( U10m - U_ice ) !
	874	! ------------------------------------------------------------ !
	875	! C-grid ice dynamics : U & V-points (same as ocean)
	876	DO jj = 2, jpjm1
	877	DO ji = fs_2, fs_jpim1 ! vect. opt.
[12015]	878	putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) &
	879	& + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) &
[11304]	880	& * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) )
[12015]	881	pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) &
	882	& + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) &
[11304]	883	& * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) )
	884	END DO
[6723]	885	END DO
[11304]	886	CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
	887	!
	888	IF(ln_ctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
	889	& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
[11334]	890	ELSE
	891	zztmp1 = 11637800.0_wp
[11360]	892	zztmp2 = -5897.8_wp
	893	DO jj = 1, jpj
[11334]	894	DO ji = 1, jpi
	895	pcd_dui(ji,jj) = zcd_dui (ji,jj)
[11360]	896	pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
	897	pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
	898	zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??)
[12015]	899	pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj)
[11334]	900	END DO
	901	END DO
	902	ENDIF
[9019]	903	!
[11304]	904	IF(ln_ctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
[9767]	905	!
[11304]	906	END SUBROUTINE blk_ice_1
[6723]	907
	908
[11304]	909	SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow )
[6723]	910	!!---------------------------------------------------------------------
[11304]	911	!! * ROUTINE blk_ice_2 *
[6723]	912	!!
[9019]	913	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
[6723]	914	!!
	915	!! ** Method : compute heat and freshwater exchanged
	916	!! between atmosphere and sea-ice using bulk formulation
	917	!! formulea, ice variables and read atmmospheric fields.
	918	!!
	919	!! caution : the net upward water flux has with mm/day unit
	920	!!---------------------------------------------------------------------
[11360]	921	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K]
[9019]	922	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
	923	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
[6727]	924	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
[11304]	925	REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair
	926	REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi
	927	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp
	928	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw
	929	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec
	930	REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
[6723]	931	!!
[6727]	932	INTEGER :: ji, jj, jl ! dummy loop indices
[9454]	933	REAL(wp) :: zst3 ! local variable
[6727]	934	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
[12015]	935	REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - -
[9019]	936	REAL(wp) :: zfr1, zfr2 ! local variables
[9454]	937	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
[9019]	938	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
	939	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
	940	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
	941	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
[9656]	942	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
[12015]	943	REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB
[6723]	944	!!---------------------------------------------------------------------
	945	!
[12015]	946	zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars
	947	zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
[6723]	948	!
[12015]	949	SELECT CASE( nhumi )
	950	CASE( np_humi_sph )
	951	zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity!
	952	CASE( np_humi_dpt )
	953	zqair(:,:) = q_sat( phumi(:,:), pslp )
	954	CASE( np_humi_rlh )
	955	zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file
	956	END SELECT
[6723]	957	!
	958	zztmp = 1. / ( 1. - albo )
[12015]	959	WHERE( ptsu(:,:,:) /= 0._wp )
	960	z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
	961	ELSEWHERE
	962	z1_st(:,:,:) = 0._wp
[9454]	963	END WHERE
[7753]	964	! ! ========================== !
	965	DO jl = 1, jpl ! Loop over ice categories !
	966	! ! ========================== !
[6723]	967	DO jj = 1 , jpj
	968	DO ji = 1, jpi
	969	! ----------------------------!
	970	! I Radiative FLUXES !
	971	! ----------------------------!
[9454]	972	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
[6723]	973	! Short Wave (sw)
	974	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
	975	! Long Wave (lw)
[12015]	976	z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
[6723]	977	! lw sensitivity
	978	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
	979
	980	! ----------------------------!
	981	! II Turbulent FLUXES !
	982	! ----------------------------!
	983
[11360]	984	! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
[6723]	985	! Sensible Heat
[12015]	986	z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))
[6723]	987	! Latent Heat
[11304]	988	zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) )
[12015]	989	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * &
	990	& ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) )
[6723]	991	! Latent heat sensitivity for ice (Dqla/Dt)
	992	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
[11360]	993	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * &
[11304]	994	& z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2
[6723]	995	ELSE
	996	dqla_ice(ji,jj,jl) = 0._wp
	997	ENDIF
	998
	999	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
[12015]	1000	z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj)
	1001
[6723]	1002	! ----------------------------!
	1003	! III Total FLUXES !
	1004	! ----------------------------!
	1005	! Downward Non Solar flux
	1006	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
	1007	! Total non solar heat flux sensitivity for ice
	1008	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
	1009	END DO
	1010	!
	1011	END DO
	1012	!
	1013	END DO
	1014	!
[11304]	1015	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
	1016	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
	1017	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
	1018	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	1019
	1020	! --- evaporation --- !
[9935]	1021	z1_rLsub = 1._wp / rLsub
	1022	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
	1023	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
	1024	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
[6723]	1025
[7753]	1026	! --- evaporation minus precipitation --- !
	1027	zsnw(:,:) = 0._wp
[9019]	1028	CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
	1029	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]	1030	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
	1031	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]	1032
[7753]	1033	! --- heat flux associated with emp --- !
[9019]	1034	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
[11304]	1035	& + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair
[7753]	1036	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
[11304]	1037	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]	1038	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
[11304]	1039	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	1040
[7753]	1041	! --- total solar and non solar fluxes --- !
[9019]	1042	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
	1043	& + qemp_ice(:,:) + qemp_oce(:,:)
	1044	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]	1045
[7753]	1046	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[11304]	1047	qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
	1048
[7504]	1049	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
	1050	DO jl = 1, jpl
[9935]	1051	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
[12015]	1052	! ! But we do not have Tice => consider it at 0degC => evap=0
[7504]	1053	END DO
	1054
[9019]	1055	! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
	1056	zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm
	1057	zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1
[6723]	1058	!
[12015]	1059	WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
[9910]	1060	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
[9019]	1061	ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
[9910]	1062	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
[9019]	1063	ELSEWHERE ! zero when hs>0
[12015]	1064	qtr_ice_top(:,:,:) = 0._wp
[9019]	1065	END WHERE
[6723]	1066	!
	1067	IF(ln_ctl) THEN
	1068	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
	1069	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
	1070	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
	1071	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
	1072	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
	1073	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
[12015]	1074	ENDIF
[6723]	1075	!
[11304]	1076	END SUBROUTINE blk_ice_2
[6723]	1077
[12015]	1078
[10531]	1079	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
[9019]	1080	!!---------------------------------------------------------------------
	1081	!! * ROUTINE blk_ice_qcn *
[6723]	1082	!!
[9019]	1083	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
	1084	!! to force sea ice / snow thermodynamics
[10534]	1085	!! in the case conduction flux is emulated
[12015]	1086	!!
[9019]	1087	!! ** Method : compute surface energy balance assuming neglecting heat storage
	1088	!! following the 0-layer Semtner (1976) approach
[6727]	1089	!!
[9019]	1090	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
	1091	!! - qcn_ice : surface inner conduction flux (W/m2)
	1092	!!
	1093	!!---------------------------------------------------------------------
[10531]	1094	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
[9076]	1095	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
	1096	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
	1097	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
	1098	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
[6723]	1099	!
[9019]	1100	INTEGER , PARAMETER :: nit = 10 ! number of iterations
	1101	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
[6723]	1102	!
[9019]	1103	INTEGER :: ji, jj, jl ! dummy loop indices
	1104	INTEGER :: iter ! local integer
	1105	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
	1106	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
	1107	REAL(wp) :: zqc, zqnet !
	1108	REAL(wp) :: zhe, zqa0 !
	1109	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
	1110	!!---------------------------------------------------------------------
[12015]	1111
[9019]	1112	! -------------------------------------!
	1113	! I Enhanced conduction factor !
	1114	! -------------------------------------!
[10531]	1115	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
[9019]	1116	! Fichefet and Morales Maqueda, JGR 1997
[6723]	1117	!
[9019]	1118	zgfac(:,:,:) = 1._wp
[12015]	1119
[10531]	1120	IF( ld_virtual_itd ) THEN
[9019]	1121	!
[9935]	1122	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
[9019]	1123	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
	1124	zfac3 = 2._wp / zepsilon
[12015]	1125	!
	1126	DO jl = 1, jpl
[9019]	1127	DO jj = 1 , jpj
	1128	DO ji = 1, jpi
[9935]	1129	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
[9019]	1130	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
	1131	END DO
	1132	END DO
	1133	END DO
[12015]	1134	!
[10531]	1135	ENDIF
[12015]	1136
[9019]	1137	! -------------------------------------------------------------!
	1138	! II Surface temperature and conduction flux !
	1139	! -------------------------------------------------------------!
[6723]	1140	!
[9935]	1141	zfac = rcnd_i * rn_cnd_s
[6723]	1142	!
[9019]	1143	DO jl = 1, jpl
	1144	DO jj = 1 , jpj
	1145	DO ji = 1, jpi
[12015]	1146	!
[9019]	1147	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
[9935]	1148	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
[9019]	1149	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
	1150	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
[9910]	1151	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
[6727]	1152	!
[9019]	1153	DO iter = 1, nit ! --- Iterative loop
	1154	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
	1155	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
	1156	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
	1157	END DO
	1158	!
	1159	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
	1160	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
	1161	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
[9910]	1162	qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) &
[12015]	1163	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]	1164
[9938]	1165	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
[12015]	1166	hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)
[9938]	1167
[9019]	1168	END DO
[6723]	1169	END DO
[9019]	1170	!
[12015]	1171	END DO
	1172	!
[9019]	1173	END SUBROUTINE blk_ice_qcn
[6723]	1174
[12015]	1175
[11360]	1176	SUBROUTINE Cdn10_Lupkes2012( pcd )
[7355]	1177	!!----------------------------------------------------------------------
	1178	!! * ROUTINE Cdn10_Lupkes2012 *
	1179	!!
[12015]	1180	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
[9019]	1181	!! to make it dependent on edges at leads, melt ponds and flows.
[7355]	1182	!! After some approximations, this can be resumed to a dependency
	1183	!! on ice concentration.
[12015]	1184	!!
[7355]	1185	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
	1186	!! with the highest level of approximation: level4, eq.(59)
	1187	!! The generic drag over a cell partly covered by ice can be re-written as follows:
	1188	!!
	1189	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
	1190	!!
	1191	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
	1192	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
	1193	!! A is the concentration of ice minus melt ponds (if any)
	1194	!!
	1195	!! This new drag has a parabolic shape (as a function of A) starting at
[12015]	1196	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
[7355]	1197	!! and going down to Cdi(say 1.4e-3) for A=1
	1198	!!
[7507]	1199	!! It is theoretically applicable to all ice conditions (not only MIZ)
[7355]	1200	!! => see Lupkes et al (2013)
	1201	!!
	1202	!! ** References : Lupkes et al. JGR 2012 (theory)
	1203	!! Lupkes et al. GRL 2013 (application to GCM)
	1204	!!
	1205	!!----------------------------------------------------------------------
[11360]	1206	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd
[7355]	1207	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
	1208	REAL(wp), PARAMETER :: znu = 1._wp
	1209	REAL(wp), PARAMETER :: zmu = 1._wp
	1210	REAL(wp), PARAMETER :: zbeta = 1._wp
	1211	REAL(wp) :: zcoef
	1212	!!----------------------------------------------------------------------
	1213	zcoef = znu + 1._wp / ( 10._wp * zbeta )
	1214
	1215	! generic drag over a cell partly covered by ice
[7507]	1216	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
	1217	!! & Cd_ice * at_i_b(:,:) + & ! pure ice drag
	1218	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
[7355]	1219
	1220	! ice-atm drag
[12015]	1221	pcd(:,:) = rCd_ice + & ! pure ice drag
[11360]	1222	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
[12015]	1223
[7355]	1224	END SUBROUTINE Cdn10_Lupkes2012
[9019]	1225
	1226
[12015]	1227	SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )
[9019]	1228	!!----------------------------------------------------------------------
	1229	!! * ROUTINE Cdn10_Lupkes2015 *
	1230	!!
	1231	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
[12015]	1232	!! between sea-ice and atmosphere with distinct momentum
	1233	!! and heat coefficients depending on sea-ice concentration
[9019]	1234	!! and atmospheric stability (no meltponds effect for now).
[12015]	1235	!!
[9019]	1236	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
	1237	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
	1238	!! it considers specific skin and form drags (Andreas et al. 2010)
[12015]	1239	!! to compute neutral transfert coefficients for both heat and
[9019]	1240	!! momemtum fluxes. Atmospheric stability effect on transfert
	1241	!! coefficient is also taken into account following Louis (1979).
	1242	!!
	1243	!! ** References : Lupkes et al. JGR 2015 (theory)
	1244	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
	1245	!!
	1246	!!----------------------------------------------------------------------
	1247	!
[11360]	1248	REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K]
[12015]	1249	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa]
[11360]	1250	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient
	1251	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient
	1252	REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat
[9019]	1253	!
	1254	! ECHAM6 constants
	1255	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
	1256	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
	1257	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
	1258	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
	1259	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
	1260	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
	1261	REAL(wp), PARAMETER :: zc2 = zc * zc
	1262	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
	1263	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
	1264	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
	1265	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
	1266	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
	1267	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
	1268	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
	1269	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
	1270	!
	1271	INTEGER :: ji, jj ! dummy loop indices
	1272	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
	1273	REAL(wp) :: zrib_o, zrib_i
	1274	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
	1275	REAL(wp) :: zChn_skin_ice, zChn_form_ice
	1276	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
	1277	REAL(wp) :: zCdn_form_tmp
[12015]	1278	!!----------------------------------------------------------------------
	1279
[9019]	1280	! Momentum Neutral Transfert Coefficients (should be a constant)
	1281	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
	1282	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
[11360]	1283	zCdn_ice = zCdn_skin_ice ! Eq. 7
[9019]	1284	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
	1285
	1286	! Heat Neutral Transfert Coefficients
[11360]	1287	zChn_skin_ice = vkarmn*2 / ( LOG( rn_zu / z0_ice + 1._wp ) LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52
[12015]	1288
[9019]	1289	! Atmospheric and Surface Variables
[10511]	1290	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
[12015]	1291	zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg]
	1292	zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg]
[9019]	1293	!
[12015]	1294	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
[9019]	1295	DO ji = fs_2, fs_jpim1
	1296	! Virtual potential temperature [K]
[10511]	1297	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
[11360]	1298	zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
[10511]	1299	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
[12015]	1300
[9019]	1301	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
	1302	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
	1303	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
[12015]	1304
[9019]	1305	! Momentum and Heat Neutral Transfert Coefficients
	1306	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
[12015]	1307	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
	1308
	1309	! Momentum and Heat Stability functions (!!GS: possibility to use psi_m_ecmwf instead ?)
[9019]	1310	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
[11360]	1311	z0i = z0_skin_ice ! over ice
[9019]	1312	IF( zrib_o <= 0._wp ) THEN
	1313	zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10
	1314	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
	1315	& )zgamma )z1_gamma
	1316	ELSE
	1317	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
	1318	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
	1319	ENDIF
[12015]	1320
[9019]	1321	IF( zrib_i <= 0._wp ) THEN
	1322	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
	1323	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
	1324	ELSE
	1325	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
	1326	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
	1327	ENDIF
[12015]	1328
[9019]	1329	! Momentum Transfert Coefficients (Eq. 38)
[11360]	1330	pcd(ji,jj) = zCdn_skin_ice * zfmi + &
	1331	& zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
[12015]	1332
[9019]	1333	! Heat Transfert Coefficients (Eq. 49)
[11360]	1334	pch(ji,jj) = zChn_skin_ice * zfhi + &
	1335	& zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
[9019]	1336	!
	1337	END DO
	1338	END DO
	1339	!
[11360]	1340	CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1., pch, 'T', 1. )
	1341	!
[9019]	1342	END SUBROUTINE Cdn10_Lupkes2015
	1343
[7355]	1344	#endif
	1345
[6723]	1346	!!======================================================================
	1347	END MODULE sbcblk

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

Download in other formats: