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_r12072_MERGE_OPTION2_2019/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk.F90 @ 12183

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

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

Download in other formats: