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.

icesbc.F90 in NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90 @ 12588

Last change on this file since 12588 was 12588, checked in by gsamson, 4 years ago

revised ABL model version including:

albmod cleaning
new abl mixing length option (nn_amxl = 3)
use rho_air function from aerobulk everywhere
remove mxl_abl (replaced by master (mxlm_abl) and dissipative (mxlm_abl) mixing lengths)
temporary flag "ln_tpot" to disable potential temperature computation in sbcblk

Property svn:keywords set to Id

File size: 17.7 KB

Rev	Line
[10535]	1	MODULE icesbc
[8586]	2	!!======================================================================
[10535]	3	!! * MODULE icesbc *
	4	!! Sea-Ice : air-ice sbc fields
[8586]	5	!!=====================================================================
[9604]	6	!! History : 4.0 ! 2017-08 (C. Rousset) Original code
	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!----------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' : SI3 sea-ice model
[8586]	12	!!----------------------------------------------------------------------
	13	USE oce ! ocean dynamics and tracers
	14	USE dom_oce ! ocean space and time domain
	15	USE ice ! sea-ice: variables
	16	USE sbc_oce ! Surface boundary condition: ocean fields
	17	USE sbc_ice ! Surface boundary condition: ice fields
	18	USE usrdef_sbc ! Surface boundary condition: user defined
	19	USE sbcblk ! Surface boundary condition: bulk
	20	USE sbccpl ! Surface boundary condition: coupled interface
[8884]	21	USE icealb ! sea-ice: albedo
[8586]	22	!
	23	USE in_out_manager ! I/O manager
	24	USE iom ! I/O manager library
	25	USE lib_mpp ! MPP library
	26	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	27	USE lbclnk ! lateral boundary conditions (or mpp links)
	28	USE timing ! Timing
[12377]	29	USE fldread !!GS: needed by agrif
[8586]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[10535]	34	PUBLIC ice_sbc_tau ! called by icestp.F90
	35	PUBLIC ice_sbc_flx ! called by icestp.F90
	36	PUBLIC ice_sbc_init ! called by icestp.F90
[8586]	37
	38	!! * Substitutions
[12377]	39	# include "do_loop_substitute.h90"
[8586]	40	!!----------------------------------------------------------------------
[10068]	41	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	42	!! $Id$
[10068]	43	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	44	!!----------------------------------------------------------------------
	45	CONTAINS
	46
[10535]	47	SUBROUTINE ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
[8586]	48	!!-------------------------------------------------------------------
[10535]	49	!! * ROUTINE ice_sbc_tau *
[8586]	50	!!
	51	!! ** Purpose : provide surface boundary condition for sea ice (momentum)
	52	!!
	53	!! ** Action : It provides the following fields:
	54	!! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
	55	!!-------------------------------------------------------------------
[9169]	56	INTEGER , INTENT(in ) :: kt ! ocean time step
	57	INTEGER , INTENT(in ) :: ksbc ! type of sbc flux
	58	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: utau_ice, vtau_ice ! air-ice stress [N/m2]
[8586]	59	!!
	60	INTEGER :: ji, jj ! dummy loop index
	61	REAL(wp), DIMENSION(jpi,jpj) :: zutau_ice, zvtau_ice
	62	!!-------------------------------------------------------------------
[9169]	63	!
[10535]	64	IF( ln_timing ) CALL timing_start('ice_sbc')
[9169]	65	!
[8586]	66	IF( kt == nit000 .AND. lwp ) THEN
	67	WRITE(numout,*)
[10535]	68	WRITE(numout,*)'ice_sbc_tau: Surface boundary condition for sea ice (momentum)'
[8586]	69	WRITE(numout,*)'~~~~~~~~~~~~~~~'
	70	ENDIF
[9169]	71	!
[8586]	72	SELECT CASE( ksbc )
	73	CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation
[12377]	74	CASE( jp_blk ) ; CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), &
	75	& sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), &
	76	& sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su , & ! inputs
	77	& putaui = utau_ice, pvtaui = vtau_ice ) ! outputs
	78	! CASE( jp_abl ) utau_ice & vtau_ice are computed in ablmod
[8920]	79	CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation
[8586]	80	END SELECT
[9169]	81	!
[8586]	82	IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation
	83	CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
[12377]	84	DO_2D_00_00
	85	utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
	86	vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
	87	END_2D
[10535]	88	CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
[8586]	89	ENDIF
[9169]	90	!
[10535]	91	IF( ln_timing ) CALL timing_stop('ice_sbc')
[8586]	92	!
[10535]	93	END SUBROUTINE ice_sbc_tau
[8586]	94
	95
[10535]	96	SUBROUTINE ice_sbc_flx( kt, ksbc )
[8586]	97	!!-------------------------------------------------------------------
[10535]	98	!! * ROUTINE ice_sbc_flx *
[8586]	99	!!
	100	!! ** Purpose : provide surface boundary condition for sea ice (flux)
	101	!!
	102	!! ** Action : It provides the following fields used in sea ice model:
	103	!! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s]
	104	!! sprecip = solid precipitation [Kg/m2/s]
	105	!! evap_ice = sublimation [Kg/m2/s]
	106	!! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2]
	107	!! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2]
	108	!! dqns_ice = non solar heat sensistivity [W/m2]
	109	!! qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
	110	!! + some fields that are not used outside this module:
	111	!! qla_ice = latent heat flux over ice [W/m2]
	112	!! dqla_ice = latent heat sensistivity [W/m2]
	113	!! tprecip = total precipitation [Kg/m2/s]
	114	!! alb_ice = albedo above sea ice
	115	!!-------------------------------------------------------------------
	116	INTEGER, INTENT(in) :: kt ! ocean time step
	117	INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk or Pure Coupled)
	118	!
[11536]	119	INTEGER :: ji, jj, jl ! dummy loop index
	120	REAL(wp) :: zmiss_val ! missing value retrieved from xios
	121	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky
	122	REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zalb, zmsk00 ! 2D workspace
[8586]	123	!!--------------------------------------------------------------------
	124	!
[10535]	125	IF( ln_timing ) CALL timing_start('ice_sbc_flx')
[8586]	126
	127	IF( kt == nit000 .AND. lwp ) THEN
	128	WRITE(numout,*)
[10535]	129	WRITE(numout,*)'ice_sbc_flx: Surface boundary condition for sea ice (flux)'
[8586]	130	WRITE(numout,*)'~~~~~~~~~~~~~~~'
	131	ENDIF
	132
[11536]	133	! get missing value from xml
	134	CALL iom_miss_val( "icetemp", zmiss_val )
	135
[8586]	136	! --- cloud-sky and overcast-sky ice albedos --- !
[8637]	137	CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
[8586]	138
	139	! albedo depends on cloud fraction because of non-linear spectral effects
	140	!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
	141	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
	142	!
[8920]	143	!
	144	SELECT CASE( ksbc ) !== fluxes over sea ice ==!
	145	!
	146	CASE( jp_usr ) !--- user defined formulation
[8962]	147	CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
[12377]	148	CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation
	149	CALL blk_ice_2 ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), &
[12588]	150	& sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )
[8906]	151	IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
[8920]	152	IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
[10534]	153	! ! compute conduction flux and surface temperature (as in Jules surface module)
	154	IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
	155	& CALL blk_ice_qcn ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
[8920]	156	CASE ( jp_purecpl ) !--- coupled formulation
[8906]	157	CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
[8920]	158	IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
[8586]	159	END SELECT
	160
[8920]	161	!--- output ice albedo and surface albedo ---!
[11536]	162	IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN
	163
	164	ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) )
	165
[11575]	166	WHERE( at_i_b < 1.e-03 )
[11536]	167	zmsk00(:,:) = 0._wp
	168	zalb (:,:) = rn_alb_oce
	169	ELSEWHERE
	170	zmsk00(:,:) = 1._wp
	171	zalb (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
[8586]	172	END WHERE
[11536]	173	! ice albedo
	174	CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )
	175	! ice+ocean albedo
[8586]	176	zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
[11536]	177	CALL iom_put( 'albedo' , zalb )
	178
	179	DEALLOCATE( zalb, zmsk00 )
	180
[8586]	181	ENDIF
	182	!
[10535]	183	IF( ln_timing ) CALL timing_stop('ice_sbc_flx')
[8586]	184	!
[10535]	185	END SUBROUTINE ice_sbc_flx
[8586]	186
	187
[10534]	188	SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
[8586]	189	!!-------------------------------------------------------------------
	190	!! * ROUTINE ice_flx_dist *
	191	!!
	192	!! ** Purpose : update the ice surface boundary condition by averaging
	193	!! and/or redistributing fluxes on ice categories
	194	!!
	195	!! ** Method : average then redistribute
	196	!!
[10534]	197	!! ** Action : depends on k_flxdist
[8586]	198	!! = -1 Do nothing (needs N(cat) fluxes)
	199	!! = 0 Average N(cat) fluxes then apply the average over the N(cat) ice
	200	!! = 1 Average N(cat) fluxes then redistribute over the N(cat) ice
	201	!! using T-ice and albedo sensitivity
	202	!! = 2 Redistribute a single flux over categories
	203	!!-------------------------------------------------------------------
[10534]	204	INTEGER , INTENT(in ) :: k_flxdist ! redistributor
[8586]	205	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
	206	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
	207	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
	208	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
	209	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
	210	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
	211	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
	212	!
	213	INTEGER :: jl ! dummy loop index
	214	!
	215	REAL(wp), DIMENSION(jpi,jpj) :: z1_at_i ! inverse of concentration
	216	!
	217	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
	218	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
	219	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories
	220	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories
	221	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
	222	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories
	223	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
	224	!!----------------------------------------------------------------------
	225	!
	226	WHERE ( at_i (:,:) > 0._wp ) ; z1_at_i(:,:) = 1._wp / at_i (:,:)
	227	ELSEWHERE ; z1_at_i(:,:) = 0._wp
	228	END WHERE
	229
[10534]	230	SELECT CASE( k_flxdist ) !== averaged on all ice categories ==!
[8586]	231	!
	232	CASE( 0 , 1 )
	233	!
	234	ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )
	235	!
	236	z_qns_m (:,:) = SUM( a_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
	237	z_qsr_m (:,:) = SUM( a_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
	238	z_dqn_m (:,:) = SUM( a_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
	239	z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
	240	z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
	241	DO jl = 1, jpl
	242	pqns_ice (:,:,jl) = z_qns_m (:,:)
	243	pqsr_ice (:,:,jl) = z_qsr_m (:,:)
	244	pdqn_ice (:,:,jl) = z_dqn_m (:,:)
	245	pevap_ice (:,:,jl) = z_evap_m(:,:)
	246	pdevap_ice(:,:,jl) = z_devap_m(:,:)
	247	END DO
	248	!
	249	DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m )
	250	!
	251	END SELECT
	252	!
[10534]	253	SELECT CASE( k_flxdist ) !== redistribution on all ice categories ==!
[8586]	254	!
	255	CASE( 1 , 2 )
	256	!
	257	ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )
	258	!
	259	zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
	260	ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
	261	DO jl = 1, jpl
	262	pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	263	pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	264	pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
	265	END DO
	266	!
	267	DEALLOCATE( zalb_m, ztem_m )
	268	!
	269	END SELECT
	270	!
	271	END SUBROUTINE ice_flx_dist
	272
[9169]	273
[10535]	274	SUBROUTINE ice_sbc_init
[8586]	275	!!-------------------------------------------------------------------
[10535]	276	!! * ROUTINE ice_sbc_init *
[8586]	277	!!
[9169]	278	!! ** Purpose : Physical constants and parameters linked to the ice dynamics
	279	!!
[10535]	280	!! ** Method : Read the namsbc namelist and check the ice-dynamic
[9169]	281	!! parameter values called at the first timestep (nit000)
[8586]	282	!!
[10535]	283	!! ** input : Namelist namsbc
[8586]	284	!!-------------------------------------------------------------------
[9169]	285	INTEGER :: ios, ioptio ! Local integer
[8586]	286	!!
[10535]	287	NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
[8586]	288	!!-------------------------------------------------------------------
	289	!
[10535]	290	READ ( numnam_ice_ref, namsbc, IOSTAT = ios, ERR = 901)
[11536]	291	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist' )
[10535]	292	READ ( numnam_ice_cfg, namsbc, IOSTAT = ios, ERR = 902 )
[11536]	293	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist' )
[10535]	294	IF(lwm) WRITE( numoni, namsbc )
[8586]	295	!
	296	IF(lwp) THEN ! control print
	297	WRITE(numout,*)
[10535]	298	WRITE(numout,*) 'ice_sbc_init: ice parameters for ice dynamics '
[9169]	299	WRITE(numout,*) '~~~~~~~~~~~~~~~~'
[10535]	300	WRITE(numout,*) ' Namelist namsbc:'
[10534]	301	WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio
	302	WRITE(numout,*) ' coefficient for ice-lead partition of snowfall rn_blow_s = ', rn_blow_s
	303	WRITE(numout,*) ' Multicategory heat flux formulation nn_flxdist = ', nn_flxdist
	304	WRITE(numout,*) ' Use conduction flux as surface condition ln_cndflx = ', ln_cndflx
	305	WRITE(numout,*) ' emulate conduction flux ln_cndemulate = ', ln_cndemulate
[8586]	306	ENDIF
	307	!
	308	IF(lwp) WRITE(numout,*)
[9570]	309	SELECT CASE( nn_flxdist ) ! SI3 Multi-category heat flux formulation
[8586]	310	CASE( -1 )
[9570]	311	IF(lwp) WRITE(numout,*) ' SI3: use per-category fluxes (nn_flxdist = -1) '
[8586]	312	CASE( 0 )
[9570]	313	IF(lwp) WRITE(numout,*) ' SI3: use average per-category fluxes (nn_flxdist = 0) '
[8586]	314	CASE( 1 )
[9570]	315	IF(lwp) WRITE(numout,*) ' SI3: use average then redistribute per-category fluxes (nn_flxdist = 1) '
	316	IF( ln_cpl ) CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in coupled mode must be /=1' )
[8586]	317	CASE( 2 )
[9570]	318	IF(lwp) WRITE(numout,*) ' SI3: Redistribute a single flux over categories (nn_flxdist = 2) '
	319	IF( .NOT. ln_cpl ) CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in forced mode must be /=2' )
[8586]	320	CASE DEFAULT
[9570]	321	CALL ctl_stop( 'ice_thd_init: SI3 option, nn_flxdist, should be between -1 and 2' )
[8586]	322	END SELECT
	323	!
[10535]	324	END SUBROUTINE ice_sbc_init
[8586]	325
	326	#else
	327	!!----------------------------------------------------------------------
[9570]	328	!! Default option : Empty module NO SI3 sea-ice model
[8586]	329	!!----------------------------------------------------------------------
	330	#endif
	331
	332	!!======================================================================
[10535]	333	END MODULE icesbc

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

Download in other formats: