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_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 12038

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

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

Download in other formats: