Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sbcblk.F90 @ 14735

Last change on this file since 14735 was 14075, checked in by cguiavarch, 4 years ago
UKMO/NEMO_4.0.4_mirror : Remove SVN keywords.
File size: 73.3 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
[9019]	17	!! ! ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
[7163]	18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
[10534]	19	!! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle)
[6723]	20	!!----------------------------------------------------------------------
	21
	22	!!----------------------------------------------------------------------
[7163]	23	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
	24	!! sbc_blk : bulk formulation as ocean surface boundary condition
[6727]	25	!! blk_oce : computes momentum, heat and freshwater fluxes over ocean
	26	!! rho_air : density of (moist) air (depends on T_air, q_air and SLP
	27	!! cp_air : specific heat of (moist) air (depends spec. hum. q_air)
	28	!! q_sat : saturation humidity as a function of SLP and temperature
	29	!! L_vap : latent heat of vaporization of water as a function of temperature
[9019]	30	!! sea-ice case only :
	31	!! blk_ice_tau : provide the air-ice stress
	32	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
[10534]	33	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	34	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
	35	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[6723]	36	!!----------------------------------------------------------------------
	37	USE oce ! ocean dynamics and tracers
	38	USE dom_oce ! ocean space and time domain
	39	USE phycst ! physical constants
	40	USE fldread ! read input fields
	41	USE sbc_oce ! Surface boundary condition: ocean fields
	42	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
	43	USE sbcdcy ! surface boundary condition: diurnal cycle
	44	USE sbcwave , ONLY : cdn_wave ! wave module
	45	USE sbc_ice ! Surface boundary condition: ice fields
	46	USE lib_fortran ! to use key_nosignedzero
[9570]	47	#if defined key_si3
[13284]	48	USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice
	49	USE icevar ! for CALL ice_var_snwblow
[6723]	50	#endif
	51	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)
	52	USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)
	53	USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
	54	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)
[6727]	55	!
[6723]	56	USE iom ! I/O manager library
	57	USE in_out_manager ! I/O manager
	58	USE lib_mpp ! distribued memory computing library
	59	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	60	USE prtctl ! Print control
	61
	62	IMPLICIT NONE
	63	PRIVATE
	64
[7163]	65	PUBLIC sbc_blk_init ! called in sbcmod
	66	PUBLIC sbc_blk ! called in sbcmod
[9570]	67	#if defined key_si3
[10535]	68	PUBLIC blk_ice_tau ! routine called in icesbc
	69	PUBLIC blk_ice_flx ! routine called in icesbc
	70	PUBLIC blk_ice_qcn ! routine called in icesbc
[9019]	71	#endif
[6723]	72
[6727]	73	!!Lolo: should ultimately be moved in the module with all physical constants ?
	74	!!gm : In principle, yes.
[6723]	75	REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg]
	76	REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg]
	77	REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
	78	REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
	79	REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
[12901]	80	REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
[6723]	81
[13284]	82	INTEGER , PARAMETER :: jpfld =11 ! maximum number of files to read
[6723]	83	INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
	84	INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
[7163]	85	INTEGER , PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
	86	INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
	87	INTEGER , PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
	88	INTEGER , PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
[6723]	89	INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
	90	INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
[7421]	91	INTEGER , PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa)
[13284]	92	INTEGER , PARAMETER :: jp_cc =10 ! index of cloud cover (-) range:0-1
	93	INTEGER , PARAMETER :: jp_tdif =11 ! index of tau diff associated to HF tau (N/m2) at T-point
[6723]	94
	95	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	96
	97	! !!! Bulk parameters
[7355]	98	REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...)
	99	REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation
	100	REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant
[9019]	101	REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice
[7355]	102	REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant
[6723]	103	!
	104	! !!* Namelist namsbc_blk : bulk parameters
	105	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
	106	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
	107	LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013)
	108	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31)
	109	!
	110	LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data
	111	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
[9019]	112	REAL(wp) :: rn_efac ! multiplication factor for evaporation
	113	REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress
[6723]	114	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
	115	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
[9019]	116	!!gm ref namelist initialize it so remove the setting to false below
	117	LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
	118	LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
[7355]	119	!
[9019]	120	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau)
	121	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens)
	122	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat)
	123	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
	124	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme)
	125	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
[6723]	126
	127	INTEGER :: nblk ! choice of the bulk algorithm
	128	! ! associated indices:
	129	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
	130	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
	131	INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)
	132	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31)
	133
	134	!! * Substitutions
	135	# include "vectopt_loop_substitute.h90"
	136	!!----------------------------------------------------------------------
[9598]	137	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	138	!! $Id$
[10068]	139	!! Software governed by the CeCILL license (see ./LICENSE)
[6723]	140	!!----------------------------------------------------------------------
	141	CONTAINS
	142
[7355]	143	INTEGER FUNCTION sbc_blk_alloc()
	144	!!-------------------------------------------------------------------
	145	!! * ROUTINE sbc_blk_alloc *
	146	!!-------------------------------------------------------------------
[9019]	147	ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
	148	& cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
[7355]	149	!
[10425]	150	CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
	151	IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
[7355]	152	END FUNCTION sbc_blk_alloc
	153
[9019]	154
[7163]	155	SUBROUTINE sbc_blk_init
	156	!!---------------------------------------------------------------------
	157	!! * ROUTINE sbc_blk_init *
	158	!!
	159	!! ** Purpose : choose and initialize a bulk formulae formulation
	160	!!
	161	!! ** Method :
	162	!!
	163	!!----------------------------------------------------------------------
[13284]	164	INTEGER :: jfpr, jfld ! dummy loop indice and argument
[7163]	165	INTEGER :: ios, ierror, ioptio ! Local integer
	166	!!
	167	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
	168	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	169	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
	170	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
[13284]	171	TYPE(FLD_N) :: sn_slp , sn_tdif, sn_cc ! " "
[7163]	172	NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields
[13284]	173	& sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, sn_cc, &
[7163]	174	& ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm
[7355]	175	& cn_dir , ln_taudif, rn_zqt, rn_zu, &
[9019]	176	& rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
[7163]	177	!!---------------------------------------------------------------------
	178	!
[7355]	179	! ! allocate sbc_blk_core array
	180	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
	181	!
[7163]	182	! !** read bulk namelist
	183	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
	184	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
[11536]	185	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
[7163]	186	!
	187	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
	188	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
[11536]	189	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
[7163]	190	!
	191	IF(lwm) WRITE( numond, namsbc_blk )
	192	!
	193	! !** initialization of the chosen bulk formulae (+ check)
	194	! !* select the bulk chosen in the namelist and check the choice
[8666]	195	ioptio = 0
[7163]	196	IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF
	197	IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF
	198	IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF
	199	IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF
	200	!
	201	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
	202	!
	203	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
[11536]	204	IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
[7163]	205	IF( sn_qsr%ln_tint ) THEN
	206	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
	207	& ' ==> We force time interpolation = .false. for qsr' )
	208	sn_qsr%ln_tint = .false.
	209	ENDIF
	210	ENDIF
	211	! !* set the bulk structure
	212	! !- store namelist information in an array
	213	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
	214	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
	215	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
	216	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
[13284]	217	slf_i(jp_slp) = sn_slp ; slf_i(jp_cc) = sn_cc
	218	slf_i(jp_tdif) = sn_tdif
[7163]	219	!
	220	lhftau = ln_taudif !- add an extra field if HF stress is used
	221	jfld = jpfld - COUNT( (/.NOT.lhftau/) )
	222	!
	223	! !- allocate the bulk structure
	224	ALLOCATE( sf(jfld), STAT=ierror )
	225	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
[8530]	226
[7163]	227	! !- fill the bulk structure with namelist informations
	228	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
	229	!
[13284]	230	DO jfpr = 1, jfld
	231	!
	232	IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero)
	233	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
	234	sf(jfpr)%fnow(:,:,1) = 0._wp
	235	ELSE !-- used field --!
	236	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
	237	IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) )
	238	IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) &
	239	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
	240	& ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
	241	ENDIF
	242	ENDDO
	243	! fill cloud cover array with constant value if "not used"
	244	IF( TRIM(sf(jp_cc)%clrootname) == 'NOT USED' ) sf(jp_cc)%fnow(:,:,1) = pp_cldf
	245
[7431]	246	IF ( ln_wave ) THEN
	247	!Activated wave module but neither drag nor stokes drift activated
[9033]	248	IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN
[10425]	249	CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
[7431]	250	!drag coefficient read from wave model definable only with mfs bulk formulae and core
	251	ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN
[10190]	252	CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
[7431]	253	ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN
	254	CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
	255	ENDIF
	256	ELSE
[9033]	257	IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &
[7431]	258	& CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &
	259	& 'with drag coefficient (ln_cdgw =T) ' , &
	260	& 'or Stokes Drift (ln_sdw=T) ' , &
[9033]	261	& 'or ocean stress modification due to waves (ln_tauwoc=T) ', &
[7431]	262	& 'or Stokes-Coriolis term (ln_stcori=T)' )
	263	ENDIF
	264	!
[7163]	265	!
	266	IF(lwp) THEN !** Control print
	267	!
	268	WRITE(numout,) ! namelist
	269	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
	270	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
	271	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
[13348]	272	WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p5
[7163]	273	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF
	274	WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif
	275	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
	276	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
	277	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
	278	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
	279	WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
	280	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
[9019]	281	WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12
	282	WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15
[7163]	283	!
	284	WRITE(numout,*)
	285	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
[9190]	286	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
	287	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
	288	CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)'
	289	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)'
[7163]	290	END SELECT
	291	!
	292	ENDIF
	293	!
	294	END SUBROUTINE sbc_blk_init
	295
	296
[6723]	297	SUBROUTINE sbc_blk( kt )
	298	!!---------------------------------------------------------------------
	299	!! * ROUTINE sbc_blk *
	300	!!
	301	!! ** Purpose : provide at each time step the surface ocean fluxes
[9019]	302	!! (momentum, heat, freshwater and runoff)
[6723]	303	!!
	304	!! ** Method : (1) READ each fluxes in NetCDF files:
	305	!! the 10m wind velocity (i-component) (m/s) at T-point
	306	!! the 10m wind velocity (j-component) (m/s) at T-point
	307	!! the 10m or 2m specific humidity ( % )
	308	!! the solar heat (W/m2)
	309	!! the Long wave (W/m2)
	310	!! the 10m or 2m air temperature (Kelvin)
	311	!! the total precipitation (rain+snow) (Kg/m2/s)
	312	!! the snow (solid prcipitation) (kg/m2/s)
	313	!! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T)
	314	!! (2) CALL blk_oce
	315	!!
	316	!! C A U T I O N : never mask the surface stress fields
	317	!! the stress is assumed to be in the (i,j) mesh referential
	318	!!
	319	!! ** Action : defined at each time-step at the air-sea interface
	320	!! - utau, vtau i- and j-component of the wind stress
	321	!! - taum wind stress module at T-point
	322	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
	323	!! - qns, qsr non-solar and solar heat fluxes
	324	!! - emp upward mass flux (evapo. - precip.)
	325	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
	326	!!
	327	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
	328	!! Brodeau et al. Ocean Modelling 2010
	329	!!----------------------------------------------------------------------
	330	INTEGER, INTENT(in) :: kt ! ocean time step
	331	!!---------------------------------------------------------------------
	332	!
	333	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
[7163]	334	!
[6723]	335	! ! compute the surface ocean fluxes using bulk formulea
	336	IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
	337
	338	#if defined key_cice
	339	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[7753]	340	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
	341	IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
	342	ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
[7282]	343	ENDIF
[7753]	344	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
	345	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
	346	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
	347	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
	348	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
	349	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
[6723]	350	ENDIF
	351	#endif
	352	!
	353	END SUBROUTINE sbc_blk
	354
	355
	356	SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
	357	!!---------------------------------------------------------------------
	358	!! * ROUTINE blk_oce *
	359	!!
	360	!! ** Purpose : provide the momentum, heat and freshwater fluxes at
	361	!! the ocean surface at each time step
	362	!!
	363	!! ** Method : bulk formulea for the ocean using atmospheric
	364	!! fields read in sbc_read
	365	!!
	366	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
	367	!! - vtau : j-component of the stress at V-point (N/m2)
	368	!! - taum : Wind stress module at T-point (N/m2)
	369	!! - wndm : Wind speed module at T-point (m/s)
	370	!! - qsr : Solar heat flux over the ocean (W/m2)
	371	!! - qns : Non Solar heat flux over the ocean (W/m2)
	372	!! - emp : evaporation minus precipitation (kg/m2/s)
	373	!!
	374	!! ** Nota : sf has to be a dummy argument for AGRIF on NEC
	375	!!---------------------------------------------------------------------
	376	INTEGER , INTENT(in ) :: kt ! time step index
	377	TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data
	378	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius]
	379	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
	380	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
	381	!
	382	INTEGER :: ji, jj ! dummy loop indices
	383	REAL(wp) :: zztmp ! local variable
[9019]	384	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
	385	REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst
	386	REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes
	387	REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation
	388	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
	389	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
	390	REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]
	391	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3]
[6723]	392	!!---------------------------------------------------------------------
	393	!
[7753]	394	! local scalars ( place there for vector optimisation purposes)
	395	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
[6723]	396
[13284]	397	! --- cloud cover --- !
	398	cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
	399
[6723]	400	! ----------------------------------------------------------------------------- !
	401	! 0 Wind components and module at T-point relative to the moving ocean !
	402	! ----------------------------------------------------------------------------- !
	403
[7753]	404	! ... components ( U10m - U_oce ) at T-point (unmasked)
	405	!!gm move zwnd_i (_j) set to zero inside the key_cyclone ???
	406	zwnd_i(:,:) = 0._wp
	407	zwnd_j(:,:) = 0._wp
[6723]	408	#if defined key_cyclone
	409	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
	410	DO jj = 2, jpjm1
	411	DO ji = fs_2, fs_jpim1 ! vect. opt.
	412	sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
	413	sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
	414	END DO
	415	END DO
	416	#endif
	417	DO jj = 2, jpjm1
	418	DO ji = fs_2, fs_jpim1 ! vect. opt.
	419	zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
	420	zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
	421	END DO
	422	END DO
[10425]	423	CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
[6723]	424	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
[7753]	425	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
	426	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
[6723]	427
	428	! ----------------------------------------------------------------------------- !
	429	! I Radiative FLUXES !
	430	! ----------------------------------------------------------------------------- !
	431
	432	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
	433	zztmp = 1. - albo
	434	IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
[7753]	435	ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
[6723]	436	ENDIF
	437
[7753]	438	zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)zst(:,:)zst(:,:)zst(:,:) ) tmask(:,:,1) ! Long Wave
[6723]	439
	440	! ----------------------------------------------------------------------------- !
	441	! II Turbulent FLUXES !
	442	! ----------------------------------------------------------------------------- !
	443
	444	! ... specific humidity at SST and IST tmask(
[6727]	445	zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
[6723]	446	!!
	447	!! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
	448	!! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
	449	!! (since reanalysis products provide T at z, not theta !)
[6727]	450	ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
[6723]	451
	452	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
[6727]	453	!
	454	CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2
[9019]	455	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]	456	CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0
[9019]	457	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]	458	CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5
[9019]	459	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]	460	CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF
[9019]	461	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6723]	462	CASE DEFAULT
[7163]	463	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
[6723]	464	END SELECT
	465
[6727]	466	! ! Compute true air density :
	467	IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrhogravrn_zu from SLP...)
[9019]	468	zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) )
[6727]	469	ELSE ! At zt:
	470	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[6723]	471	END IF
	472
[9019]	473	!! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef.
	474	!! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef.
[7355]	475
[6727]	476	DO jj = 1, jpj ! tau module, i and j component
[6723]	477	DO ji = 1, jpi
[9019]	478	zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed
[6723]	479	taum (ji,jj) = zztmp * wndm (ji,jj)
	480	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
	481	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
	482	END DO
	483	END DO
	484
[6727]	485	! ! add the HF tau contribution to the wind stress module
[7753]	486	IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
[6727]	487
[6723]	488	CALL iom_put( "taum_oce", taum ) ! output wind stress module
	489
	490	! ... utau, vtau at U- and V_points, resp.
	491	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
	492	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
	493	DO jj = 1, jpjm1
	494	DO ji = 1, fs_jpim1
	495	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
	496	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
	497	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
	498	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
	499	END DO
	500	END DO
[10425]	501	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
[6723]	502
	503	! Turbulent fluxes over ocean
	504	! -----------------------------
	505
	506	! zqla used as temporary array, for rho*U (common term of bulk formulae):
[9727]	507	zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
[6723]	508
	509	IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
	510	!! q_air and t_air are given at 10m (wind reference height)
[9019]	511	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
	512	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed
[6723]	513	ELSE
	514	!! q_air and t_air are not given at 10m (wind reference height)
	515	! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
[9019]	516	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
	517	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed
[6723]	518	ENDIF
	519
[6727]	520	zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux
[6723]	521
	522
	523	IF(ln_ctl) THEN
[9019]	524	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' )
	525	CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' )
[6723]	526	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
	527	CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' )
	528	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
	529	& tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask )
	530	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ')
	531	CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ')
	532	ENDIF
	533
	534	! ----------------------------------------------------------------------------- !
	535	! III Total FLUXES !
	536	! ----------------------------------------------------------------------------- !
	537	!
[7753]	538	emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.)
	539	& - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)
	540	!
	541	qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar
[9935]	542	& - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
[7753]	543	& - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST
	544	& + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair
	545	& * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &
	546	& + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
[9935]	547	& * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
[9727]	548	qns(:,:) = qns(:,:) * tmask(:,:,1)
[7753]	549	!
[9570]	550	#if defined key_si3
[9656]	551	qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by SI3)
[7753]	552	qsr_oce(:,:) = qsr(:,:)
[6723]	553	#endif
	554	!
	555	IF ( nn_ice == 0 ) THEN
	556	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
	557	CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean
	558	CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean
	559	CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean
	560	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
	561	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
	562	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
[9727]	563	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
	564	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
[9019]	565	CALL iom_put( 'snowpre', sprecip ) ! Snow
	566	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	567	ENDIF
	568	!
	569	IF(ln_ctl) THEN
	570	CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ')
	571	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
	572	CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ')
	573	CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
	574	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
	575	ENDIF
	576	!
	577	END SUBROUTINE blk_oce
	578
	579
[9019]	580
	581	FUNCTION rho_air( ptak, pqa, pslp )
	582	!!-------------------------------------------------------------------------------
	583	!! * FUNCTION rho_air *
	584	!!
	585	!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
	586	!!
	587	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
	588	!!-------------------------------------------------------------------------------
	589	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	590	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
	591	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]
	592	REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]
	593	!!-------------------------------------------------------------------------------
	594	!
	595	rho_air = pslp / ( R_dryptak ( 1._wp + rctv0*pqa ) )
	596	!
	597	END FUNCTION rho_air
	598
	599
	600	FUNCTION cp_air( pqa )
	601	!!-------------------------------------------------------------------------------
	602	!! * FUNCTION cp_air *
	603	!!
	604	!! ** Purpose : Compute specific heat (Cp) of moist air
	605	!!
	606	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
	607	!!-------------------------------------------------------------------------------
	608	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
	609	REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg]
	610	!!-------------------------------------------------------------------------------
	611	!
	612	Cp_air = Cp_dry + Cp_vap * pqa
	613	!
	614	END FUNCTION cp_air
	615
	616
	617	FUNCTION q_sat( ptak, pslp )
	618	!!----------------------------------------------------------------------------------
	619	!! * FUNCTION q_sat *
	620	!!
	621	!! ** Purpose : Specific humidity at saturation in [kg/kg]
	622	!! Based on accurate estimate of "e_sat"
	623	!! aka saturation water vapor (Goff, 1957)
	624	!!
	625	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
	626	!!----------------------------------------------------------------------------------
	627	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	628	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
	629	REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]
	630	!
	631	INTEGER :: ji, jj ! dummy loop indices
	632	REAL(wp) :: ze_sat, ztmp ! local scalar
	633	!!----------------------------------------------------------------------------------
	634	!
	635	DO jj = 1, jpj
	636	DO ji = 1, jpi
	637	!
	638	ztmp = rt0 / ptak(ji,jj)
	639	!
	640	! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
	641	ze_sat = 10.*( 10.79574(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &
	642	& + 1.5047510.(-4)(1. - 10.*(-8.2969(ptak(ji,jj)/rt0 - 1.)) ) &
	643	& + 0.4287310.(-3)(10.*(4.76955(1. - ztmp)) - 1.) + 0.78614 )
	644	!
	645	q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wppslp(ji,jj) - (1._wp - reps0)ze_sat ) ! 0.01 because SLP is in [Pa]
	646	!
	647	END DO
	648	END DO
	649	!
	650	END FUNCTION q_sat
	651
	652
	653	FUNCTION gamma_moist( ptak, pqa )
	654	!!----------------------------------------------------------------------------------
	655	!! * FUNCTION gamma_moist *
	656	!!
	657	!! ** Purpose : Compute the moist adiabatic lapse-rate.
	658	!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
	659	!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
	660	!!
	661	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
	662	!!----------------------------------------------------------------------------------
	663	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	664	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]
	665	REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate
	666	!
	667	INTEGER :: ji, jj ! dummy loop indices
	668	REAL(wp) :: zrv, ziRT ! local scalar
	669	!!----------------------------------------------------------------------------------
	670	!
	671	DO jj = 1, jpj
	672	DO ji = 1, jpi
	673	zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
	674	ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT
[9935]	675	gamma_moist(ji,jj) = grav * ( 1. + rLevapzrvziRT ) / ( Cp_dry + rLevaprLevapzrvreps0ziRT/ptak(ji,jj) )
[9019]	676	END DO
	677	END DO
	678	!
	679	END FUNCTION gamma_moist
	680
	681
	682	FUNCTION L_vap( psst )
	683	!!---------------------------------------------------------------------------------
	684	!! * FUNCTION L_vap *
	685	!!
	686	!! ** Purpose : Compute the latent heat of vaporization of water from temperature
	687	!!
	688	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
	689	!!----------------------------------------------------------------------------------
	690	REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]
	691	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
	692	!!----------------------------------------------------------------------------------
	693	!
	694	L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6
	695	!
	696	END FUNCTION L_vap
	697
[9570]	698	#if defined key_si3
[9019]	699	!!----------------------------------------------------------------------
[9570]	700	!! 'key_si3' SI3 sea-ice model
[9019]	701	!!----------------------------------------------------------------------
	702	!! blk_ice_tau : provide the air-ice stress
	703	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
[10534]	704	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	705	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
	706	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
	707	!!----------------------------------------------------------------------
	708
[6723]	709	SUBROUTINE blk_ice_tau
	710	!!---------------------------------------------------------------------
	711	!! * ROUTINE blk_ice_tau *
	712	!!
	713	!! ** Purpose : provide the surface boundary condition over sea-ice
	714	!!
	715	!! ** Method : compute momentum using bulk formulation
	716	!! formulea, ice variables and read atmospheric fields.
	717	!! NB: ice drag coefficient is assumed to be a constant
	718	!!---------------------------------------------------------------------
	719	INTEGER :: ji, jj ! dummy loop indices
[9019]	720	REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point
	721	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
[12926]	722	REAL(wp) :: zztmp1 , zztmp2 ! temporary values
[9019]	723	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau)
[6723]	724	!!---------------------------------------------------------------------
	725	!
[9019]	726	! set transfer coefficients to default sea-ice values
	727	Cd_atm(:,:) = Cd_ice
	728	Ch_atm(:,:) = Cd_ice
	729	Ce_atm(:,:) = Cd_ice
[6723]	730
[9019]	731	wndm_ice(:,:) = 0._wp !!gm brutal....
[7355]	732
[9019]	733	! ------------------------------------------------------------ !
	734	! Wind module relative to the moving ice ( U10m - U_ice ) !
	735	! ------------------------------------------------------------ !
[9767]	736	! C-grid ice dynamics : U & V-points (same as ocean)
	737	DO jj = 2, jpjm1
	738	DO ji = fs_2, fs_jpim1 ! vect. opt.
	739	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
	740	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
	741	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
[9019]	742	END DO
[9767]	743	END DO
[10425]	744	CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
[9767]	745	!
[9019]	746	! Make ice-atm. drag dependent on ice concentration
	747	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
	748	CALL Cdn10_Lupkes2012( Cd_atm )
	749	Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical
	750	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
	751	CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
[7355]	752	ENDIF
	753
[9019]	754	!! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.
	755	!! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef.
	756
[6723]	757	! local scalars ( place there for vector optimisation purposes)
	758	! Computing density of air! Way denser that 1.2 over sea-ice !!!
[7355]	759	zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
[6723]	760
[9019]	761	! ------------------------------------------------------------ !
	762	! Wind stress relative to the moving ice ( U10m - U_ice ) !
	763	! ------------------------------------------------------------ !
[12926]	764	zztmp1 = rn_vfac * 0.5_wp
	765	DO jj = 2, jpj ! at T point
	766	DO ji = 2, jpi
	767	zztmp2 = zrhoa(ji,jj) * Cd_atm(ji,jj) * wndm_ice(ji,jj)
	768	utau_ice(ji,jj) = zztmp2 * ( sf(jp_wndi)%fnow(ji,jj,1) - zztmp1 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
	769	vtau_ice(ji,jj) = zztmp2 * ( sf(jp_wndj)%fnow(ji,jj,1) - zztmp1 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
	770	END DO
	771	END DO
	772	!
	773	DO jj = 2, jpjm1 ! U & V-points (same as ocean).
[9767]	774	DO ji = fs_2, fs_jpim1 ! vect. opt.
[12926]	775	! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
	776	zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
	777	zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
	778	utau_ice(ji,jj) = zztmp1 * ( utau_ice(ji,jj) + utau_ice(ji+1,jj ) )
	779	vtau_ice(ji,jj) = zztmp2 * ( vtau_ice(ji,jj) + vtau_ice(ji ,jj+1) )
[6723]	780	END DO
[9767]	781	END DO
[10425]	782	CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
[9019]	783	!
[9767]	784	!
[6723]	785	IF(ln_ctl) THEN
	786	CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
	787	CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
	788	ENDIF
[9019]	789	!
[6723]	790	END SUBROUTINE blk_ice_tau
	791
	792
[9019]	793	SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
[6723]	794	!!---------------------------------------------------------------------
	795	!! * ROUTINE blk_ice_flx *
	796	!!
[9019]	797	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
[6723]	798	!!
	799	!! ** Method : compute heat and freshwater exchanged
	800	!! between atmosphere and sea-ice using bulk formulation
	801	!! formulea, ice variables and read atmmospheric fields.
	802	!!
	803	!! caution : the net upward water flux has with mm/day unit
	804	!!---------------------------------------------------------------------
[6727]	805	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature
[9019]	806	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
	807	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
[6727]	808	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
[6723]	809	!!
[6727]	810	INTEGER :: ji, jj, jl ! dummy loop indices
[9454]	811	REAL(wp) :: zst3 ! local variable
[6727]	812	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
[13284]	813	REAL(wp) :: zztmp, z1_rLsub ! - -
[9454]	814	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
[9019]	815	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
	816	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
	817	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
	818	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
[9656]	819	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
[9019]	820	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa
[12276]	821	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
[13284]	822	REAL(wp), DIMENSION(jpi,jpj) :: ztri
[6723]	823	!!---------------------------------------------------------------------
	824	!
[9019]	825	zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars
	826	zcoef_dqla = -Ls * 11637800. * (-5897.8)
[6723]	827	!
[6727]	828	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[6723]	829	!
	830	zztmp = 1. / ( 1. - albo )
[9454]	831	WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
	832	ELSEWHERE ; z1_st(:,:,:) = 0._wp
	833	END WHERE
[7753]	834	! ! ========================== !
	835	DO jl = 1, jpl ! Loop over ice categories !
	836	! ! ========================== !
[6723]	837	DO jj = 1 , jpj
	838	DO ji = 1, jpi
	839	! ----------------------------!
	840	! I Radiative FLUXES !
	841	! ----------------------------!
[9454]	842	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
[6723]	843	! Short Wave (sw)
	844	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
	845	! Long Wave (lw)
	846	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
	847	! lw sensitivity
	848	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
	849
	850	! ----------------------------!
	851	! II Turbulent FLUXES !
	852	! ----------------------------!
	853
[9019]	854	! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
[6723]	855	! Sensible Heat
[9019]	856	z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
[6723]	857	! Latent Heat
[9019]	858	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
[9454]	859	& ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
[6723]	860	! Latent heat sensitivity for ice (Dqla/Dt)
	861	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
[9454]	862	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
	863	& z1_st(ji,jj,jl)z1_st(ji,jj,jl) EXP(-5897.8 * z1_st(ji,jj,jl))
[6723]	864	ELSE
	865	dqla_ice(ji,jj,jl) = 0._wp
	866	ENDIF
	867
	868	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
[9019]	869	z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
[6723]	870
	871	! ----------------------------!
	872	! III Total FLUXES !
	873	! ----------------------------!
	874	! Downward Non Solar flux
	875	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
	876	! Total non solar heat flux sensitivity for ice
	877	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
	878	END DO
	879	!
	880	END DO
	881	!
	882	END DO
	883	!
[9727]	884	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
	885	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
[9019]	886	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
	887	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	888
	889	! --- evaporation --- !
[9935]	890	z1_rLsub = 1._wp / rLsub
	891	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
	892	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
	893	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
[6723]	894
[7753]	895	! --- evaporation minus precipitation --- !
	896	zsnw(:,:) = 0._wp
[13284]	897	CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
[9019]	898	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]	899	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
	900	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]	901
[7753]	902	! --- heat flux associated with emp --- !
[9019]	903	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
[7753]	904	& + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair
	905	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
[9935]	906	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]	907	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
[9935]	908	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	909
[7753]	910	! --- total solar and non solar fluxes --- !
[9019]	911	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
	912	& + qemp_ice(:,:) + qemp_oce(:,:)
	913	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]	914
[7753]	915	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[9935]	916	qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	917
[7504]	918	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
	919	DO jl = 1, jpl
[9935]	920	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
[9019]	921	! ! But we do not have Tice => consider it at 0degC => evap=0
[7504]	922	END DO
	923
[13284]	924	! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
	925	IF( nn_qtrice == 0 ) THEN
	926	! formulation derived from Grenfell and Maykut (1977), where transmission rate
	927	! 1) depends on cloudiness
	928	! 2) is 0 when there is any snow
	929	! 3) tends to 1 for thin ice
	930	ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm
	931	DO jl = 1, jpl
	932	WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
	933	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
	934	ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
	935	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
	936	ELSEWHERE ! zero when hs>0
	937	qtr_ice_top(:,:,jl) = 0._wp
	938	END WHERE
	939	ENDDO
	940	ELSEIF( nn_qtrice == 1 ) THEN
	941	! formulation is derived from the thesis of M. Lebrun (2019).
	942	! It represents the best fit using several sets of observations
	943	! It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
	944	qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
	945	ENDIF
[6723]	946	!
[12276]	947
	948	IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
	949	ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
	950	CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
	951	CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average)
	952	ENDIF
	953	IF( iom_use('hflx_rain_cea') ) THEN
	954	ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
	955	CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average)
	956	ENDIF
	957	IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN
	958	WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ; ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
	959	ELSEWHERE ; ztmp(:,:) = rcp * sst_m(:,:)
	960	ENDWHERE
	961	ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
	962	CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average)
	963	CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
	964	CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice)
	965	ENDIF
	966	!
[6723]	967	IF(ln_ctl) THEN
	968	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
	969	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
	970	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
	971	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
	972	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
	973	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
	974	ENDIF
	975	!
	976	END SUBROUTINE blk_ice_flx
[6727]	977
[6723]	978
[10531]	979	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
[9019]	980	!!---------------------------------------------------------------------
	981	!! * ROUTINE blk_ice_qcn *
[6723]	982	!!
[9019]	983	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
	984	!! to force sea ice / snow thermodynamics
[10534]	985	!! in the case conduction flux is emulated
[9019]	986	!!
	987	!! ** Method : compute surface energy balance assuming neglecting heat storage
	988	!! following the 0-layer Semtner (1976) approach
[6727]	989	!!
[9019]	990	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
	991	!! - qcn_ice : surface inner conduction flux (W/m2)
	992	!!
	993	!!---------------------------------------------------------------------
[10531]	994	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
[9076]	995	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
	996	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
	997	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
	998	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
[6723]	999	!
[9019]	1000	INTEGER , PARAMETER :: nit = 10 ! number of iterations
	1001	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
[6723]	1002	!
[9019]	1003	INTEGER :: ji, jj, jl ! dummy loop indices
	1004	INTEGER :: iter ! local integer
	1005	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
	1006	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
	1007	REAL(wp) :: zqc, zqnet !
	1008	REAL(wp) :: zhe, zqa0 !
	1009	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
	1010	!!---------------------------------------------------------------------
	1011
	1012	! -------------------------------------!
	1013	! I Enhanced conduction factor !
	1014	! -------------------------------------!
[10531]	1015	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
[9019]	1016	! Fichefet and Morales Maqueda, JGR 1997
[6723]	1017	!
[9019]	1018	zgfac(:,:,:) = 1._wp
	1019
[10531]	1020	IF( ld_virtual_itd ) THEN
[9019]	1021	!
[9935]	1022	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
[9019]	1023	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
	1024	zfac3 = 2._wp / zepsilon
	1025	!
	1026	DO jl = 1, jpl
	1027	DO jj = 1 , jpj
	1028	DO ji = 1, jpi
[9935]	1029	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
[9019]	1030	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
	1031	END DO
	1032	END DO
	1033	END DO
	1034	!
[10531]	1035	ENDIF
[9019]	1036
	1037	! -------------------------------------------------------------!
	1038	! II Surface temperature and conduction flux !
	1039	! -------------------------------------------------------------!
[6723]	1040	!
[9935]	1041	zfac = rcnd_i * rn_cnd_s
[6723]	1042	!
[9019]	1043	DO jl = 1, jpl
	1044	DO jj = 1 , jpj
	1045	DO ji = 1, jpi
	1046	!
	1047	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
[9935]	1048	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
[9019]	1049	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
	1050	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
[9910]	1051	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
[6727]	1052	!
[9019]	1053	DO iter = 1, nit ! --- Iterative loop
	1054	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
	1055	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
	1056	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
	1057	END DO
	1058	!
	1059	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
	1060	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
	1061	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
[9910]	1062	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) ) &
[9019]	1063	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]	1064
[9938]	1065	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
	1066	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)
	1067
[9019]	1068	END DO
[6723]	1069	END DO
[9019]	1070	!
	1071	END DO
	1072	!
	1073	END SUBROUTINE blk_ice_qcn
	1074
[6723]	1075
[7355]	1076	SUBROUTINE Cdn10_Lupkes2012( Cd )
	1077	!!----------------------------------------------------------------------
	1078	!! * ROUTINE Cdn10_Lupkes2012 *
	1079	!!
[9019]	1080	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
	1081	!! to make it dependent on edges at leads, melt ponds and flows.
[7355]	1082	!! After some approximations, this can be resumed to a dependency
	1083	!! on ice concentration.
	1084	!!
	1085	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
	1086	!! with the highest level of approximation: level4, eq.(59)
	1087	!! The generic drag over a cell partly covered by ice can be re-written as follows:
	1088	!!
	1089	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
	1090	!!
	1091	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
	1092	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
	1093	!! A is the concentration of ice minus melt ponds (if any)
	1094	!!
	1095	!! This new drag has a parabolic shape (as a function of A) starting at
	1096	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
	1097	!! and going down to Cdi(say 1.4e-3) for A=1
	1098	!!
[7507]	1099	!! It is theoretically applicable to all ice conditions (not only MIZ)
[7355]	1100	!! => see Lupkes et al (2013)
	1101	!!
	1102	!! ** References : Lupkes et al. JGR 2012 (theory)
	1103	!! Lupkes et al. GRL 2013 (application to GCM)
	1104	!!
	1105	!!----------------------------------------------------------------------
	1106	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
	1107	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
	1108	REAL(wp), PARAMETER :: znu = 1._wp
	1109	REAL(wp), PARAMETER :: zmu = 1._wp
	1110	REAL(wp), PARAMETER :: zbeta = 1._wp
	1111	REAL(wp) :: zcoef
	1112	!!----------------------------------------------------------------------
	1113	zcoef = znu + 1._wp / ( 10._wp * zbeta )
	1114
	1115	! generic drag over a cell partly covered by ice
[7507]	1116	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
	1117	!! & Cd_ice * at_i_b(:,:) + & ! pure ice drag
	1118	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
[7355]	1119
	1120	! ice-atm drag
[7507]	1121	Cd(:,:) = Cd_ice + & ! pure ice drag
	1122	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
[7355]	1123
	1124	END SUBROUTINE Cdn10_Lupkes2012
[9019]	1125
	1126
	1127	SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
	1128	!!----------------------------------------------------------------------
	1129	!! * ROUTINE Cdn10_Lupkes2015 *
	1130	!!
	1131	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
	1132	!! between sea-ice and atmosphere with distinct momentum
	1133	!! and heat coefficients depending on sea-ice concentration
	1134	!! and atmospheric stability (no meltponds effect for now).
	1135	!!
	1136	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
	1137	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
	1138	!! it considers specific skin and form drags (Andreas et al. 2010)
	1139	!! to compute neutral transfert coefficients for both heat and
	1140	!! momemtum fluxes. Atmospheric stability effect on transfert
	1141	!! coefficient is also taken into account following Louis (1979).
	1142	!!
	1143	!! ** References : Lupkes et al. JGR 2015 (theory)
	1144	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
	1145	!!
	1146	!!----------------------------------------------------------------------
	1147	!
	1148	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
	1149	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch
[10511]	1150	REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat
[9019]	1151	!
	1152	! ECHAM6 constants
	1153	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
	1154	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
	1155	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
	1156	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
	1157	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
	1158	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
	1159	REAL(wp), PARAMETER :: zc2 = zc * zc
	1160	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
	1161	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
	1162	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
	1163	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
	1164	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
	1165	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
	1166	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
	1167	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
	1168	!
	1169	INTEGER :: ji, jj ! dummy loop indices
	1170	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
	1171	REAL(wp) :: zrib_o, zrib_i
	1172	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
	1173	REAL(wp) :: zChn_skin_ice, zChn_form_ice
	1174	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
	1175	REAL(wp) :: zCdn_form_tmp
	1176	!!----------------------------------------------------------------------
	1177
[10511]	1178	! mean temperature
	1179	WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
	1180	ELSEWHERE ; ztm_su(:,:) = rt0
	1181	ENDWHERE
	1182
[9019]	1183	! Momentum Neutral Transfert Coefficients (should be a constant)
	1184	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
	1185	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
	1186	zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)
	1187	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
	1188
	1189	! Heat Neutral Transfert Coefficients
	1190	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 (cf Lupkes email for details)
	1191
	1192	! Atmospheric and Surface Variables
[10511]	1193	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
	1194	zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg]
	1195	zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg]
[9019]	1196	!
	1197	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
	1198	DO ji = fs_2, fs_jpim1
	1199	! Virtual potential temperature [K]
[10511]	1200	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
	1201	zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
	1202	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
[9019]	1203
	1204	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
	1205	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
	1206	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
	1207
	1208	! Momentum and Heat Neutral Transfert Coefficients
	1209	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
	1210	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
	1211
	1212	! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
	1213	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
	1214	z0i = z0_skin_ice ! over ice (cf Lupkes email for details)
	1215	IF( zrib_o <= 0._wp ) THEN
	1216	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
	1217	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
	1218	& )zgamma )z1_gamma
	1219	ELSE
	1220	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
	1221	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
	1222	ENDIF
	1223
	1224	IF( zrib_i <= 0._wp ) THEN
	1225	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
	1226	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
	1227	ELSE
	1228	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
	1229	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
	1230	ENDIF
	1231
	1232	! Momentum Transfert Coefficients (Eq. 38)
	1233	Cd(ji,jj) = zCdn_skin_ice * zfmi + &
	1234	& 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) )
	1235
	1236	! Heat Transfert Coefficients (Eq. 49)
	1237	Ch(ji,jj) = zChn_skin_ice * zfhi + &
	1238	& 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) )
	1239	!
	1240	END DO
	1241	END DO
[10425]	1242	CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. )
[9019]	1243	!
	1244	END SUBROUTINE Cdn10_Lupkes2015
	1245
[7355]	1246	#endif
	1247
[6723]	1248	!!======================================================================
	1249	END MODULE sbcblk

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_mirror/src/OCE/SBC/sbcblk.F90 @ 14735

Download in other formats: