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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90 @ 12199

Last change on this file since 12199 was 12199, checked in by laurent, 4 years ago
Catches up with SBCBLK bugfixes of branch "dev_r12072_MERGE_OPTION2_2019"
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
[12182]	17	!! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
[7163]	18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
[12182]	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
[12182]	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
[12182]	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
[12182]	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
[12182]	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
[12182]	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
[12182]	68	PUBLIC blk_oce_1 ! called in sbcabl
	69	PUBLIC blk_oce_2 ! called in sbcabl
[9570]	70	#if defined key_si3
[12182]	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
[12182]	74	#endif
[6723]	75
[12182]	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
[12182]	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)
[12182]	94	LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013)
	95	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1)
[6723]	96	!
[12182]	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	!
[12182]	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
[12182]	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)
[12182]	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	!!-------------------------------------------------------------------
[12182]	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	!!
[12182]	157	!! ** Method :
[7163]	158	!!
	159	!!----------------------------------------------------------------------
[12182]	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
[12182]	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 ! " "
[12182]	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
[12182]	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	!
[12182]	179	! !** read bulk namelist
[12199]	180	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
[7163]	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	!
[12199]	184	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
[7163]	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
[12182]	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' )
[12182]	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' )
[12182]	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
[12182]	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
[12182]	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
[12182]	260	ALLOCATE( sf(jpfld), STAT=ierror )
[7163]	261	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
[12182]	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	!
[12182]	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' )
[12182]	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
[12182]	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	!
[12182]	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	!
[12182]	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
[12182]	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)'
[12182]	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	!
[12182]	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	!!
[12182]	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
[12182]	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
[12182]	400
	401	! Sanity/consistence test on humidity at first time step to detect potential screw-up:
	402	IF( kt == nit000 ) THEN
[12199]	403	IF(lwp) WRITE(numout,*) ''
[12182]	404	#if defined key_agrif
[12199]	405	IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ==='
[12182]	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
[12199]	419	IF (lwp) WRITE(numout,'(" Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp
[12182]	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
[12199]	426	IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
[12182]	427	#endif
[12199]	428	IF(lwp) WRITE(numout,*) ''
[12182]	429	END IF !IF( kt == nit000 )
[6723]	430	! ! compute the surface ocean fluxes using bulk formulea
[12182]	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)
[12199]	436	& tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out
[6723]	437
[12182]	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
[12199]	440	& sf(jp_snow)%fnow(:,:,1), tsk_m, & ! <<= in
[12182]	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)
[12182]	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)
[12182]	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
[12182]	473	SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp
[12199]	474	& pslp , pst , pu , pv, & ! inp
	475	& pqsr , pqlw , & ! inp
	476	& ptsk, pssq , pcd_du, psen , pevp ) ! out
[6723]	477	!!---------------------------------------------------------------------
[12182]	478	!! * ROUTINE blk_oce_1 *
[6723]	479	!!
[12182]	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	!!
[12182]	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	!!
[12182]	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	!!---------------------------------------------------------------------
[12182]	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]
[12199]	498	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celsius]
[12182]	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 !
[12199]	503	REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius]
[12182]	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]
[12182]	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)
[12199]	523	! ! Temporary conversion 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)
[12182]	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.
[12182]	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.
[12182]	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	! ----------------------------------------------------------------------------- !
[12182]	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
[12182]	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	! ----------------------------------------------------------------------------- !
[12182]	567	! II Turbulent FLUXES !
[6723]	568	! ----------------------------------------------------------------------------- !
	569
[12182]	570	! specific humidity at SST
[12199]	571	pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
[6723]	572
[12182]	573	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
[12199]	574	!! Backup "bulk SST" and associated spec. hum.
	575	zztmp1(:,:) = ptsk(:,:)
[12182]	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	!
[12182]	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 )
[12199]	614	CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, &
[12182]	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 )
[12199]	618	CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
[12182]	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 )
[12199]	623	CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
[12182]	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 )
[12199]	628	CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
[12182]	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' )
[12182]	634
[6723]	635	END SELECT
	636
[12182]	637	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
[12199]	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(:,:)
[12182]	645	END WHERE
[6723]	646	END IF
	647
[12182]	648	!! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef.
	649	!! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef.
[7355]	650
[12182]	651	IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
	652	!! If zu == zt, then ensuring once for all that:
	653	t_zu(:,:) = ztpot(:,:)
	654	q_zu(:,:) = zqair(:,:)
	655	ENDIF
[6723]	656
[6727]	657
[12182]	658	! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90
	659	! -------------------------------------------------------------
[6723]	660
[12182]	661	IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
	662	!! FL do we need this multiplication by tmask ... ???
	663	DO jj = 1, jpj
	664	DO ji = 1, jpi
	665	zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1)
	666	wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
	667	pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
	668	psen(ji,jj) = zztmp * zch_oce(ji,jj)
	669	pevp(ji,jj) = zztmp * zce_oce(ji,jj)
	670	END DO
[6723]	671	END DO
[12182]	672	ELSE !== BLK formulation ==! turbulent fluxes computation
[12199]	673	CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
[12182]	674	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
	675	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
	676	& taum(:,:), psen(:,:), zqla(:,:), &
	677	& pEvap=pevp(:,:), prhoa=rhoa(:,:) )
[6723]	678
[12182]	679	zqla(:,:) = zqla(:,:) * tmask(:,:,1)
	680	psen(:,:) = psen(:,:) * tmask(:,:,1)
	681	taum(:,:) = taum(:,:) * tmask(:,:,1)
	682	pevp(:,:) = pevp(:,:) * tmask(:,:,1)
[6723]	683
[12182]	684	! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array...
	685	zcd_oce = 0._wp
	686	WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm
	687	zwnd_i = zcd_oce * zwnd_i
	688	zwnd_j = zcd_oce * zwnd_j
[6723]	689
[12182]	690	CALL iom_put( "taum_oce", taum ) ! output wind stress module
[6723]	691
[12182]	692	! ... utau, vtau at U- and V_points, resp.
	693	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
	694	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
	695	DO jj = 1, jpjm1
	696	DO ji = 1, fs_jpim1
	697	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
	698	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
	699	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
	700	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
	701	END DO
	702	END DO
	703	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
[6723]	704
[12182]	705	IF(ln_ctl) THEN
	706	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
	707	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
	708	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
	709	ENDIF
	710	!
[12199]	711	ENDIF !IF( ln_abl )
	712
	713	ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
	714
	715	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
	716	CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
	717	CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
[6723]	718	ENDIF
[12199]	719
[6723]	720	IF(ln_ctl) THEN
[12182]	721	CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' )
	722	CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' )
	723	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' )
[6723]	724	ENDIF
	725	!
[12182]	726	END SUBROUTINE blk_oce_1
[6723]	727
	728
[12182]	729	SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in
[12199]	730	& psnow, ptsk, psen, pevp ) ! <<= in
[12182]	731	!!---------------------------------------------------------------------
	732	!! * ROUTINE blk_oce_2 *
[9019]	733	!!
[12182]	734	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
	735	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
	736	!! atmospheric variables (from ABL or external data)
[9019]	737	!!
[12182]	738	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
	739	!! - vtau : j-component of the stress at V-point (N/m2)
	740	!! - taum : Wind stress module at T-point (N/m2)
	741	!! - wndm : Wind speed module at T-point (m/s)
	742	!! - qsr : Solar heat flux over the ocean (W/m2)
	743	!! - qns : Non Solar heat flux over the ocean (W/m2)
	744	!! - emp : evaporation minus precipitation (kg/m2/s)
	745	!!---------------------------------------------------------------------
	746	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair
	747	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr
	748	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw
	749	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
	750	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
[12199]	751	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius]
[12182]	752	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
	753	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
[9019]	754	!
[12182]	755	INTEGER :: ji, jj ! dummy loop indices
	756	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
[12199]	757	REAL(wp), DIMENSION(jpi,jpj) :: ztskk ! skin temp. in Kelvin
	758	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes
[12182]	759	REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation
	760	!!---------------------------------------------------------------------
[9019]	761	!
[12182]	762	! local scalars ( place there for vector optimisation purposes)
[9019]	763
[12199]	764	ztskk(:,:) = ptsk(:,:) + rt0 ! => ptsk in Kelvin rather than Celsius
	765
[12182]	766	! ----------------------------------------------------------------------------- !
	767	! III Net longwave radiative FLUX !
	768	! ----------------------------------------------------------------------------- !
[9019]	769
[12182]	770	!! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
[12199]	771	!! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
	772	zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefanztskk(:,:)ztskk(:,:)ztskk(:,:)ztskk(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux
[9019]	773
[12199]	774	! Latent flux over ocean
	775	! -----------------------
[12182]	776
	777	! use scalar version of L_vap() for AGRIF compatibility
[9019]	778	DO jj = 1, jpj
	779	DO ji = 1, jpi
[12199]	780	zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
[12182]	781	ENDDO
	782	ENDDO
[9019]	783
[12182]	784	IF(ln_ctl) THEN
	785	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' )
	786	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
[9019]	787
[12182]	788	ENDIF
	789
	790	! ----------------------------------------------------------------------------- !
	791	! IV Total FLUXES !
	792	! ----------------------------------------------------------------------------- !
[9019]	793	!
[12182]	794	emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.)
	795	& - pprec(:,:) * rn_pfac ) * tmask(:,:,1)
[9019]	796	!
[12182]	797	qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar
	798	& - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
[12199]	799	& - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST
[12182]	800	& + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair
	801	& * ( ptair(:,:) - rt0 ) * rcp &
	802	& + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
	803	& * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
	804	qns(:,:) = qns(:,:) * tmask(:,:,1)
[9019]	805	!
[12182]	806	#if defined key_si3
	807	qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3)
	808	qsr_oce(:,:) = qsr(:,:)
	809	#endif
[9019]	810	!
[12182]	811	CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3]
	812	CALL iom_put( "evap_oce" , pevp ) ! evaporation
	813	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
	814	CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean
	815	CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean
	816	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
	817	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
	818	CALL iom_put( 'snowpre', sprecip ) ! Snow
	819	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[9019]	820	!
[12182]	821	IF ( nn_ice == 0 ) THEN
	822	CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean
	823	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
	824	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
	825	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
	826	ENDIF
	827	!
	828	IF(ln_ctl) THEN
	829	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
	830	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
	831	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
	832	ENDIF
	833	!
	834	END SUBROUTINE blk_oce_2
[9019]	835
[12182]	836
[9570]	837	#if defined key_si3
[9019]	838	!!----------------------------------------------------------------------
[9570]	839	!! 'key_si3' SI3 sea-ice model
[9019]	840	!!----------------------------------------------------------------------
[12182]	841	!! blk_ice_1 : provide the air-ice stress
	842	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
[10534]	843	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	844	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
[12182]	845	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[9019]	846	!!----------------------------------------------------------------------
	847
[12182]	848	SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs
	849	& putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs
[6723]	850	!!---------------------------------------------------------------------
[12182]	851	!! * ROUTINE blk_ice_1 *
[6723]	852	!!
	853	!! ** Purpose : provide the surface boundary condition over sea-ice
	854	!!
	855	!! ** Method : compute momentum using bulk formulation
	856	!! formulea, ice variables and read atmospheric fields.
	857	!! NB: ice drag coefficient is assumed to be a constant
	858	!!---------------------------------------------------------------------
[12182]	859	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
	860	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
	861	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
	862	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
	863	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s]
	864	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
	865	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
	866	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
	867	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk
	868	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk
	869	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl
	870	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl
	871	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl
	872	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
	873	!
[6723]	874	INTEGER :: ji, jj ! dummy loop indices
[9019]	875	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
[12182]	876	REAL(wp) :: zootm_su ! sea-ice surface mean temperature
	877	REAL(wp) :: zztmp1, zztmp2 ! temporary arrays
	878	REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau)
[6723]	879	!!---------------------------------------------------------------------
	880	!
	881
[9019]	882	! ------------------------------------------------------------ !
	883	! Wind module relative to the moving ice ( U10m - U_ice ) !
	884	! ------------------------------------------------------------ !
[9767]	885	! C-grid ice dynamics : U & V-points (same as ocean)
	886	DO jj = 2, jpjm1
	887	DO ji = fs_2, fs_jpim1 ! vect. opt.
[12182]	888	zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) )
	889	zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
[9767]	890	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
[9019]	891	END DO
[9767]	892	END DO
[10425]	893	CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
[9767]	894	!
[9019]	895	! Make ice-atm. drag dependent on ice concentration
	896	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
[12182]	897	CALL Cdn10_Lupkes2012( Cd_ice )
	898	Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical
	899	Ce_ice(:,:) = Cd_ice(:,:)
[9019]	900	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
[12182]	901	CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice )
	902	Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical
[7355]	903	ENDIF
	904
[12182]	905	!! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef.
	906	!! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef.
[9019]	907
[6723]	908	! local scalars ( place there for vector optimisation purposes)
[12182]	909	!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)
	910	zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
[6723]	911
[12182]	912	IF( ln_blk ) THEN
	913	! ------------------------------------------------------------ !
	914	! Wind stress relative to the moving ice ( U10m - U_ice ) !
	915	! ------------------------------------------------------------ !
	916	! C-grid ice dynamics : U & V-points (same as ocean)
	917	DO jj = 2, jpjm1
	918	DO ji = fs_2, fs_jpim1 ! vect. opt.
	919	putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) &
	920	& + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) &
	921	& * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) )
	922	pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) &
	923	& + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) &
	924	& * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) )
	925	END DO
[6723]	926	END DO
[12182]	927	CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
	928	!
	929	IF(ln_ctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
	930	& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
	931	ELSE
	932	zztmp1 = 11637800.0_wp
	933	zztmp2 = -5897.8_wp
	934	DO jj = 1, jpj
	935	DO ji = 1, jpi
	936	pcd_dui(ji,jj) = zcd_dui (ji,jj)
	937	pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
	938	pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
	939	zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??)
	940	pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj)
	941	END DO
	942	END DO
	943	ENDIF
[9019]	944	!
[12182]	945	IF(ln_ctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
[9767]	946	!
[12182]	947	END SUBROUTINE blk_ice_1
[6723]	948
	949
[12182]	950	SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow )
[6723]	951	!!---------------------------------------------------------------------
[12182]	952	!! * ROUTINE blk_ice_2 *
[6723]	953	!!
[9019]	954	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
[6723]	955	!!
	956	!! ** Method : compute heat and freshwater exchanged
	957	!! between atmosphere and sea-ice using bulk formulation
	958	!! formulea, ice variables and read atmmospheric fields.
	959	!!
	960	!! caution : the net upward water flux has with mm/day unit
	961	!!---------------------------------------------------------------------
[12182]	962	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K]
[9019]	963	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
	964	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
[6727]	965	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
[12182]	966	REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair
	967	REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi
	968	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp
	969	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw
	970	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec
	971	REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
[6723]	972	!!
[6727]	973	INTEGER :: ji, jj, jl ! dummy loop indices
[9454]	974	REAL(wp) :: zst3 ! local variable
[6727]	975	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
[12182]	976	REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - -
[9019]	977	REAL(wp) :: zfr1, zfr2 ! local variables
[9454]	978	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
[9019]	979	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
	980	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
	981	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
	982	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
[9656]	983	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
[12182]	984	REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB
[12193]	985	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
[6723]	986	!!---------------------------------------------------------------------
	987	!
[12182]	988	zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars
	989	zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
[6723]	990	!
[12182]	991	SELECT CASE( nhumi )
	992	CASE( np_humi_sph )
	993	zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity!
	994	CASE( np_humi_dpt )
	995	zqair(:,:) = q_sat( phumi(:,:), pslp )
	996	CASE( np_humi_rlh )
	997	zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file
	998	END SELECT
[6723]	999	!
	1000	zztmp = 1. / ( 1. - albo )
[12182]	1001	WHERE( ptsu(:,:,:) /= 0._wp )
	1002	z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
	1003	ELSEWHERE
	1004	z1_st(:,:,:) = 0._wp
[9454]	1005	END WHERE
[7753]	1006	! ! ========================== !
	1007	DO jl = 1, jpl ! Loop over ice categories !
	1008	! ! ========================== !
[6723]	1009	DO jj = 1 , jpj
	1010	DO ji = 1, jpi
	1011	! ----------------------------!
	1012	! I Radiative FLUXES !
	1013	! ----------------------------!
[9454]	1014	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
[6723]	1015	! Short Wave (sw)
	1016	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
	1017	! Long Wave (lw)
[12182]	1018	z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
[6723]	1019	! lw sensitivity
	1020	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
	1021
	1022	! ----------------------------!
	1023	! II Turbulent FLUXES !
	1024	! ----------------------------!
	1025
[12182]	1026	! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
[6723]	1027	! Sensible Heat
[12182]	1028	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]	1029	! Latent Heat
[12182]	1030	zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) )
	1031	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * &
	1032	& ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) )
[6723]	1033	! Latent heat sensitivity for ice (Dqla/Dt)
	1034	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
[12182]	1035	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * &
	1036	& z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2
[6723]	1037	ELSE
	1038	dqla_ice(ji,jj,jl) = 0._wp
	1039	ENDIF
	1040
	1041	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
[12182]	1042	z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj)
[6723]	1043
	1044	! ----------------------------!
	1045	! III Total FLUXES !
	1046	! ----------------------------!
	1047	! Downward Non Solar flux
	1048	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
	1049	! Total non solar heat flux sensitivity for ice
	1050	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
	1051	END DO
	1052	!
	1053	END DO
	1054	!
	1055	END DO
	1056	!
[12182]	1057	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
	1058	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
	1059	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
	1060	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	1061
	1062	! --- evaporation --- !
[9935]	1063	z1_rLsub = 1._wp / rLsub
	1064	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
	1065	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
	1066	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
[6723]	1067
[7753]	1068	! --- evaporation minus precipitation --- !
	1069	zsnw(:,:) = 0._wp
[9019]	1070	CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
	1071	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]	1072	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
	1073	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]	1074
[7753]	1075	! --- heat flux associated with emp --- !
[9019]	1076	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
[12182]	1077	& + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair
[7753]	1078	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
[12182]	1079	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]	1080	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
[12182]	1081	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	1082
[7753]	1083	! --- total solar and non solar fluxes --- !
[9019]	1084	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
	1085	& + qemp_ice(:,:) + qemp_oce(:,:)
	1086	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]	1087
[7753]	1088	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[12182]	1089	qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	1090
[7504]	1091	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
	1092	DO jl = 1, jpl
[9935]	1093	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
[12182]	1094	! ! But we do not have Tice => consider it at 0degC => evap=0
[7504]	1095	END DO
	1096
[9019]	1097	! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
	1098	zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm
	1099	zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1
[6723]	1100	!
[12182]	1101	WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
[9910]	1102	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
[9019]	1103	ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
[9910]	1104	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
[9019]	1105	ELSEWHERE ! zero when hs>0
[12182]	1106	qtr_ice_top(:,:,:) = 0._wp
[9019]	1107	END WHERE
[6723]	1108	!
[12193]	1109
	1110	IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
[12199]	1111	ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
[12193]	1112	IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
	1113	IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average)
	1114	ENDIF
	1115	IF( iom_use('hflx_rain_cea') ) THEN
	1116	ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
	1117	IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average)
	1118	ENDIF
	1119	IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN
[12199]	1120	WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 )
	1121	ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
	1122	ELSEWHERE
	1123	ztmp(:,:) = rcp * sst_m(:,:)
	1124	ENDWHERE
	1125	ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
	1126	IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average)
	1127	IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
	1128	IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice)
[12193]	1129	ENDIF
	1130	!
[6723]	1131	IF(ln_ctl) THEN
	1132	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
	1133	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
	1134	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
	1135	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
	1136	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
	1137	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
	1138	ENDIF
	1139	!
[12182]	1140	END SUBROUTINE blk_ice_2
[6723]	1141
[12182]	1142
[10531]	1143	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
[9019]	1144	!!---------------------------------------------------------------------
	1145	!! * ROUTINE blk_ice_qcn *
[6723]	1146	!!
[9019]	1147	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
	1148	!! to force sea ice / snow thermodynamics
[10534]	1149	!! in the case conduction flux is emulated
[12182]	1150	!!
[9019]	1151	!! ** Method : compute surface energy balance assuming neglecting heat storage
	1152	!! following the 0-layer Semtner (1976) approach
[6727]	1153	!!
[9019]	1154	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
	1155	!! - qcn_ice : surface inner conduction flux (W/m2)
	1156	!!
	1157	!!---------------------------------------------------------------------
[10531]	1158	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
[9076]	1159	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
	1160	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
	1161	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
	1162	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
[6723]	1163	!
[9019]	1164	INTEGER , PARAMETER :: nit = 10 ! number of iterations
	1165	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
[6723]	1166	!
[9019]	1167	INTEGER :: ji, jj, jl ! dummy loop indices
	1168	INTEGER :: iter ! local integer
	1169	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
	1170	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
	1171	REAL(wp) :: zqc, zqnet !
	1172	REAL(wp) :: zhe, zqa0 !
	1173	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
	1174	!!---------------------------------------------------------------------
[12182]	1175
[9019]	1176	! -------------------------------------!
	1177	! I Enhanced conduction factor !
	1178	! -------------------------------------!
[10531]	1179	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
[9019]	1180	! Fichefet and Morales Maqueda, JGR 1997
[6723]	1181	!
[9019]	1182	zgfac(:,:,:) = 1._wp
[12182]	1183
[10531]	1184	IF( ld_virtual_itd ) THEN
[9019]	1185	!
[9935]	1186	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
[9019]	1187	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
	1188	zfac3 = 2._wp / zepsilon
[12182]	1189	!
	1190	DO jl = 1, jpl
[9019]	1191	DO jj = 1 , jpj
	1192	DO ji = 1, jpi
[9935]	1193	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
[9019]	1194	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
	1195	END DO
	1196	END DO
	1197	END DO
[12182]	1198	!
[10531]	1199	ENDIF
[12182]	1200
[9019]	1201	! -------------------------------------------------------------!
	1202	! II Surface temperature and conduction flux !
	1203	! -------------------------------------------------------------!
[6723]	1204	!
[9935]	1205	zfac = rcnd_i * rn_cnd_s
[6723]	1206	!
[9019]	1207	DO jl = 1, jpl
	1208	DO jj = 1 , jpj
	1209	DO ji = 1, jpi
[12182]	1210	!
[9019]	1211	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
[9935]	1212	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
[9019]	1213	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
	1214	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
[9910]	1215	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
[6727]	1216	!
[9019]	1217	DO iter = 1, nit ! --- Iterative loop
	1218	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
	1219	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
	1220	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
	1221	END DO
	1222	!
	1223	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
	1224	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
	1225	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
[9910]	1226	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) ) &
[12182]	1227	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]	1228
[9938]	1229	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
[12182]	1230	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]	1231
[9019]	1232	END DO
[6723]	1233	END DO
[9019]	1234	!
[12182]	1235	END DO
	1236	!
[9019]	1237	END SUBROUTINE blk_ice_qcn
[6723]	1238
[12182]	1239
	1240	SUBROUTINE Cdn10_Lupkes2012( pcd )
[7355]	1241	!!----------------------------------------------------------------------
	1242	!! * ROUTINE Cdn10_Lupkes2012 *
	1243	!!
[12182]	1244	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
[9019]	1245	!! to make it dependent on edges at leads, melt ponds and flows.
[7355]	1246	!! After some approximations, this can be resumed to a dependency
	1247	!! on ice concentration.
[12182]	1248	!!
[7355]	1249	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
	1250	!! with the highest level of approximation: level4, eq.(59)
	1251	!! The generic drag over a cell partly covered by ice can be re-written as follows:
	1252	!!
	1253	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
	1254	!!
	1255	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
	1256	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
	1257	!! A is the concentration of ice minus melt ponds (if any)
	1258	!!
	1259	!! This new drag has a parabolic shape (as a function of A) starting at
[12182]	1260	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
[7355]	1261	!! and going down to Cdi(say 1.4e-3) for A=1
	1262	!!
[7507]	1263	!! It is theoretically applicable to all ice conditions (not only MIZ)
[7355]	1264	!! => see Lupkes et al (2013)
	1265	!!
	1266	!! ** References : Lupkes et al. JGR 2012 (theory)
	1267	!! Lupkes et al. GRL 2013 (application to GCM)
	1268	!!
	1269	!!----------------------------------------------------------------------
[12182]	1270	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd
[7355]	1271	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
	1272	REAL(wp), PARAMETER :: znu = 1._wp
	1273	REAL(wp), PARAMETER :: zmu = 1._wp
	1274	REAL(wp), PARAMETER :: zbeta = 1._wp
	1275	REAL(wp) :: zcoef
	1276	!!----------------------------------------------------------------------
	1277	zcoef = znu + 1._wp / ( 10._wp * zbeta )
	1278
	1279	! generic drag over a cell partly covered by ice
[7507]	1280	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
	1281	!! & Cd_ice * at_i_b(:,:) + & ! pure ice drag
	1282	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
[7355]	1283
	1284	! ice-atm drag
[12182]	1285	pcd(:,:) = rCd_ice + & ! pure ice drag
	1286	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
	1287
[7355]	1288	END SUBROUTINE Cdn10_Lupkes2012
[9019]	1289
	1290
[12182]	1291	SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )
[9019]	1292	!!----------------------------------------------------------------------
	1293	!! * ROUTINE Cdn10_Lupkes2015 *
	1294	!!
	1295	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
[12182]	1296	!! between sea-ice and atmosphere with distinct momentum
	1297	!! and heat coefficients depending on sea-ice concentration
[9019]	1298	!! and atmospheric stability (no meltponds effect for now).
[12182]	1299	!!
[9019]	1300	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
	1301	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
	1302	!! it considers specific skin and form drags (Andreas et al. 2010)
[12182]	1303	!! to compute neutral transfert coefficients for both heat and
[9019]	1304	!! momemtum fluxes. Atmospheric stability effect on transfert
	1305	!! coefficient is also taken into account following Louis (1979).
	1306	!!
	1307	!! ** References : Lupkes et al. JGR 2015 (theory)
	1308	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
	1309	!!
	1310	!!----------------------------------------------------------------------
	1311	!
[12182]	1312	REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K]
	1313	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa]
	1314	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient
	1315	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient
	1316	REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat
[9019]	1317	!
	1318	! ECHAM6 constants
	1319	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
	1320	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
	1321	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
	1322	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
	1323	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
	1324	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
	1325	REAL(wp), PARAMETER :: zc2 = zc * zc
	1326	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
	1327	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
	1328	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
	1329	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
	1330	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
	1331	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
	1332	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
	1333	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
	1334	!
	1335	INTEGER :: ji, jj ! dummy loop indices
	1336	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
	1337	REAL(wp) :: zrib_o, zrib_i
	1338	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
	1339	REAL(wp) :: zChn_skin_ice, zChn_form_ice
	1340	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
	1341	REAL(wp) :: zCdn_form_tmp
	1342	!!----------------------------------------------------------------------
	1343
	1344	! Momentum Neutral Transfert Coefficients (should be a constant)
	1345	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
	1346	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
[12182]	1347	zCdn_ice = zCdn_skin_ice ! Eq. 7
[9019]	1348	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
	1349
	1350	! Heat Neutral Transfert Coefficients
[12182]	1351	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
	1352
[9019]	1353	! Atmospheric and Surface Variables
[10511]	1354	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
[12182]	1355	zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg]
	1356	zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg]
[9019]	1357	!
	1358	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
	1359	DO ji = fs_2, fs_jpim1
	1360	! Virtual potential temperature [K]
[10511]	1361	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
[12182]	1362	zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
[10511]	1363	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
[12182]	1364
[9019]	1365	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
	1366	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
	1367	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
[12182]	1368
[9019]	1369	! Momentum and Heat Neutral Transfert Coefficients
	1370	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
[12182]	1371	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
	1372
	1373	! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?)
[9019]	1374	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
[12182]	1375	z0i = z0_skin_ice ! over ice
[9019]	1376	IF( zrib_o <= 0._wp ) THEN
	1377	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
	1378	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
	1379	& )zgamma )z1_gamma
	1380	ELSE
	1381	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
	1382	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
	1383	ENDIF
[12182]	1384
[9019]	1385	IF( zrib_i <= 0._wp ) THEN
	1386	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
	1387	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
	1388	ELSE
	1389	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
	1390	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
	1391	ENDIF
[12182]	1392
[9019]	1393	! Momentum Transfert Coefficients (Eq. 38)
[12182]	1394	pcd(ji,jj) = zCdn_skin_ice * zfmi + &
[9019]	1395	& 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) )
[12182]	1396
[9019]	1397	! Heat Transfert Coefficients (Eq. 49)
[12182]	1398	pch(ji,jj) = zChn_skin_ice * zfhi + &
[9019]	1399	& 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) )
	1400	!
	1401	END DO
	1402	END DO
[12182]	1403	CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1., pch, 'T', 1. )
[9019]	1404	!
	1405	END SUBROUTINE Cdn10_Lupkes2015
	1406
[7355]	1407	#endif
	1408
[6723]	1409	!!======================================================================
	1410	END MODULE sbcblk

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

Download in other formats: