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.

sbcice_cice.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8324

Last change on this file since 8324 was 8324, checked in by clem, 7 years ago
STEP3 (2): clean separation between sea-ice and ocean
Property svn:keywords set to `Id`
File size: 44.3 KB

Rev	Line
[2874]	1	MODULE sbcice_cice
	2	!!======================================================================
	3	!! * MODULE sbcice_cice *
	4	!! To couple with sea ice model CICE (LANL)
	5	!!=====================================================================
	6	#if defined key_cice
	7	!!----------------------------------------------------------------------
	8	!! 'key_cice' : CICE sea-ice model
	9	!!----------------------------------------------------------------------
	10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
	11	!!----------------------------------------------------------------------
	12	USE oce ! ocean dynamics and tracers
	13	USE dom_oce ! ocean space and time domain
[3275]	14	USE domvvl
[3625]	15	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
[2874]	16	USE in_out_manager ! I/O manager
[4990]	17	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
[2874]	18	USE lib_mpp ! distributed memory computing library
	19	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3186]	20	USE wrk_nemo ! work arrays
[3193]	21	USE timing ! Timing
[2874]	22	USE daymod ! calendar
	23	USE fldread ! read input fields
	24	USE sbc_oce ! Surface boundary condition: ocean fields
	25	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	26	USE sbcblk ! Surface boundary condition: bulk
[2874]	27	USE sbccpl
	28
	29	USE ice_kinds_mod
	30	USE ice_blocks
	31	USE ice_domain
	32	USE ice_domain_size
	33	USE ice_boundary
	34	USE ice_constants
	35	USE ice_gather_scatter
	36	USE ice_calendar, only: dt
[3625]	37	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
[4990]	38	# if defined key_cice4
[2874]	39	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	40	strocnxT,strocnyT, &
[3189]	41	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
	42	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
[2874]	43	flatn_f,fsurfn_f,fcondtopn_f, &
	44	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
	45	swvdr,swvdf,swidr,swidf
[4990]	46	USE ice_therm_vertical, only: calc_Tsfc
	47	#else
	48	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	49	strocnxT,strocnyT, &
[4990]	50	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
	51	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
	52	flatn_f,fsurfn_f,fcondtopn_f, &
	53	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
	54	swvdr,swvdf,swidr,swidf
	55	USE ice_therm_shared, only: calc_Tsfc
	56	#endif
[2874]	57	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
[3176]	58	USE ice_atmo, only: calc_strair
[2874]	59
	60	USE CICE_InitMod
	61	USE CICE_RunMod
	62	USE CICE_FinalMod
	63
	64	IMPLICIT NONE
	65	PRIVATE
	66
	67	PUBLIC cice_sbc_init ! routine called by sbc_init
	68	PUBLIC cice_sbc_final ! routine called by sbc_final
	69	PUBLIC sbc_ice_cice ! routine called by sbc
	70
[4627]	71	INTEGER :: ji_off
	72	INTEGER :: jj_off
[3625]	73
[2874]	74	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
	75	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
	76	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
	77	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
	78	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
	79	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
	80	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
	81	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
	82	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
	83	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
	84	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
	85	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
	86	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
	87	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
	88	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	89
	90	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
	91
[5836]	92	!!----------------------------------------------------------------------
	93	!! NEMO/OPA 3.7 , NEMO-consortium (2015)
[5215]	94	!! $Id$
[5836]	95	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	96	!!----------------------------------------------------------------------
[2874]	97	CONTAINS
	98
	99	INTEGER FUNCTION sbc_ice_cice_alloc()
	100	!!----------------------------------------------------------------------
	101	!! * FUNCTION sbc_ice_cice_alloc *
	102	!!----------------------------------------------------------------------
	103	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
	104	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
	105	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
	106	END FUNCTION sbc_ice_cice_alloc
	107
[4990]	108	SUBROUTINE sbc_ice_cice( kt, ksbc )
[2874]	109	!!---------------------------------------------------------------------
	110	!! * ROUTINE sbc_ice_cice *
	111	!!
	112	!! ** Purpose : update the ocean surface boundary condition via the
	113	!! CICE Sea Ice Model time stepping
	114	!!
[3040]	115	!! ** Method : - Get any extra forcing fields for CICE
	116	!! - Prepare forcing fields
[2874]	117	!! - CICE model time stepping
	118	!! - call the routine that computes mass and
	119	!! heat fluxes at the ice/ocean interface
	120	!!
	121	!! ** Action : - time evolution of the CICE sea-ice model
	122	!! - update all sbc variables below sea-ice:
[3625]	123	!! utau, vtau, qns , qsr, emp , sfx
[2874]	124	!!---------------------------------------------------------------------
	125	INTEGER, INTENT(in) :: kt ! ocean time step
[4990]	126	INTEGER, INTENT(in) :: ksbc ! surface forcing type
[2874]	127	!!----------------------------------------------------------------------
[3193]	128	!
	129	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_cice')
	130	!
[2874]	131	! !----------------------!
	132	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
	133	! !----------------------!
	134
	135	! Make sure any fluxes required for CICE are set
[4990]	136	IF ( ksbc == jp_flx ) THEN
[2874]	137	CALL cice_sbc_force(kt)
[5407]	138	ELSE IF ( ksbc == jp_purecpl ) THEN
[8316]	139	CALL sbc_cpl_ice_flx( fr_i )
[2874]	140	ENDIF
	141
[4990]	142	CALL cice_sbc_in ( kt, ksbc )
[2874]	143	CALL CICE_Run
[4990]	144	CALL cice_sbc_out ( kt, ksbc )
[2874]	145
[5407]	146	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
[2874]	147
	148	ENDIF ! End sea-ice time step only
[3193]	149	!
	150	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_cice')
[2874]	151
	152	END SUBROUTINE sbc_ice_cice
	153
[5836]	154
	155	SUBROUTINE cice_sbc_init( ksbc )
[2874]	156	!!---------------------------------------------------------------------
	157	!! * ROUTINE cice_sbc_init *
[3040]	158	!! ** Purpose: Initialise ice related fields for NEMO and coupling
[2874]	159	!!
[5836]	160	!!---------------------------------------------------------------------
[4990]	161	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[3625]	162	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
	163	REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
[4990]	164	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[2874]	165	!!---------------------------------------------------------------------
	166
[3193]	167	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_init')
	168	!
[3625]	169	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
	170	!
[2874]	171	IF(lwp) WRITE(numout,*)'cice_sbc_init'
	172
[4627]	173	ji_off = INT ( (jpiglo - nx_global) / 2 )
	174	jj_off = INT ( (jpjglo - ny_global) / 2 )
	175
[4990]	176	#if defined key_nemocice_decomp
	177	! Pass initial SST from NEMO to CICE so ice is initialised correctly if
	178	! there is no restart file.
	179	! Values from a CICE restart file would overwrite this
	180	IF ( .NOT. ln_rstart ) THEN
	181	CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)
	182	ENDIF
	183	#endif
	184
[2874]	185	! Initialize CICE
[3176]	186	CALL CICE_Initialize
[2874]	187
[3176]	188	! Do some CICE consistency checks
[5407]	189	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3193]	190	IF ( calc_strair .OR. calc_Tsfc ) THEN
	191	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
	192	ENDIF
[7646]	193	ELSEIF (ksbc == jp_blk) THEN
[3193]	194	IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
	195	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
	196	ENDIF
	197	ENDIF
[3176]	198
	199
[2874]	200	! allocate sbc_ice and sbc_cice arrays
	201	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
	202
	203	! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
	204	IF( .NOT. ln_rstart ) THEN
	205	tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
	206	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
	207	ENDIF
	208
[3193]	209	fr_iu(:,:)=0.0
	210	fr_iv(:,:)=0.0
[2874]	211
[3193]	212	CALL cice2nemo(aice,fr_i, 'T', 1. )
[5407]	213	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3625]	214	DO jl=1,ncat
	215	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	216	ENDDO
	217	ENDIF
[2874]	218
	219	! T point to U point
	220	! T point to V point
[3193]	221	DO jj=1,jpjm1
	222	DO ji=1,jpim1
	223	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	224	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	225	ENDDO
	226	ENDDO
[2874]	227
[3193]	228	CALL lbc_lnk ( fr_iu , 'U', 1. )
	229	CALL lbc_lnk ( fr_iv , 'V', 1. )
[3625]	230
[8313]	231	! set the snow+ice mass
	232	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	233	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	234	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	235	snwice_mass_b(:,:) = snwice_mass(:,:)
	236
[6140]	237	IF( .NOT.ln_rstart ) THEN
[8313]	238	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[4990]	239	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	240	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[6140]	241
	242	!!gm This should be put elsewhere.... (same remark for limsbc)
	243	!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
	244	IF( .NOT.ln_linssh ) THEN
	245	!
	246	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	247	e3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	248	e3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	249	ENDDO
	250	e3t_a(:,:,:) = e3t_b(:,:,:)
	251	! Reconstruction of all vertical scale factors at now and before time-steps
	252	! =============================================================================
	253	! Horizontal scale factor interpolations
	254	! --------------------------------------
	255	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	256	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	257	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	258	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	259	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	260	! Vertical scale factor interpolations
	261	! ------------------------------------
	262	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	263	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	264	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	265	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	266	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	267	! t- and w- points depth
	268	! ----------------------
	269	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	270	gdepw_n(:,:,1) = 0.0_wp
	271	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	272	DO jk = 2, jpk
	273	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
	274	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	275	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	276	END DO
	277	ENDIF
[4990]	278	ENDIF
[3625]	279	ENDIF
[6140]	280	!
[3625]	281	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	282	!
	283	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
	284	!
[2874]	285	END SUBROUTINE cice_sbc_init
	286
[3152]	287
[5836]	288	SUBROUTINE cice_sbc_in( kt, ksbc )
[2874]	289	!!---------------------------------------------------------------------
	290	!! * ROUTINE cice_sbc_in *
[3040]	291	!! ** Purpose: Set coupling fields and pass to CICE
[2874]	292	!!---------------------------------------------------------------------
[3152]	293	INTEGER, INTENT(in ) :: kt ! ocean time step
[4990]	294	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
[5836]	295	!
[3625]	296	INTEGER :: ji, jj, jl ! dummy loop indices
	297	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
[3152]	298	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
[3625]	299	REAL(wp) :: zintb, zintn ! dummy argument
[3152]	300	!!---------------------------------------------------------------------
[2874]	301
[3193]	302	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
	303	!
[3625]	304	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
[3152]	305	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
[2874]	306
[3193]	307	IF( kt == nit000 ) THEN
[2874]	308	IF(lwp) WRITE(numout,*)'cice_sbc_in'
[3193]	309	ENDIF
[2874]	310
[3193]	311	ztmp(:,:)=0.0
[2874]	312
	313	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
	314	! the first time-step)
	315
	316	! forced and coupled case
	317
[5407]	318	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[2874]	319
[3193]	320	ztmpn(:,:,:)=0.0
[2874]	321
	322	! x comp of wind stress (CI_1)
	323	! U point to F point
[3193]	324	DO jj=1,jpjm1
	325	DO ji=1,jpi
	326	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
	327	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
	328	ENDDO
	329	ENDDO
	330	CALL nemo2cice(ztmp,strax,'F', -1. )
[2874]	331
	332	! y comp of wind stress (CI_2)
	333	! V point to F point
[3193]	334	DO jj=1,jpj
	335	DO ji=1,jpim1
	336	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	337	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	338	ENDDO
	339	ENDDO
	340	CALL nemo2cice(ztmp,stray,'F', -1. )
[2874]	341
	342	! Surface downward latent heat flux (CI_5)
[4990]	343	IF (ksbc == jp_flx) THEN
[3625]	344	DO jl=1,ncat
	345	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
[3193]	346	ENDDO
	347	ELSE
[2874]	348	! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
[3193]	349	qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
[2874]	350	! End of temporary code
[3193]	351	DO jj=1,jpj
	352	DO ji=1,jpi
	353	IF (fr_i(ji,jj).eq.0.0) THEN
[3625]	354	DO jl=1,ncat
	355	ztmpn(ji,jj,jl)=0.0
[3193]	356	ENDDO
	357	! This will then be conserved in CICE
	358	ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
	359	ELSE
[3625]	360	DO jl=1,ncat
	361	ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
[3193]	362	ENDDO
	363	ENDIF
	364	ENDDO
	365	ENDDO
	366	ENDIF
[3625]	367	DO jl=1,ncat
	368	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
[2874]	369
	370	! GBM conductive flux through ice (CI_6)
	371	! Convert to GBM
[4990]	372	IF (ksbc == jp_flx) THEN
[3625]	373	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
[3193]	374	ELSE
[3625]	375	ztmp(:,:) = botmelt(:,:,jl)
[3193]	376	ENDIF
[3625]	377	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
[2874]	378
	379	! GBM surface heat flux (CI_7)
	380	! Convert to GBM
[4990]	381	IF (ksbc == jp_flx) THEN
[3625]	382	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
[3193]	383	ELSE
[3625]	384	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
[3193]	385	ENDIF
[3625]	386	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
[3193]	387	ENDDO
[2874]	388
[7646]	389	ELSE IF (ksbc == jp_blk) THEN
[2874]	390
[7646]	391	! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
[2874]	392	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	393	ztmp(:,:) = wndi_ice(:,:)
	394	CALL nemo2cice(ztmp,uatm,'T', -1. )
	395	ztmp(:,:) = wndj_ice(:,:)
	396	CALL nemo2cice(ztmp,vatm,'T', -1. )
	397	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
	398	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
	399	ztmp(:,:) = qsr_ice(:,:,1)
	400	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
	401	ztmp(:,:) = qlw_ice(:,:,1)
	402	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
	403	ztmp(:,:) = tatm_ice(:,:)
	404	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
	405	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
[2874]	406	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	407	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	408	! Constant (101000.) atm pressure assumed
	409	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
	410	ztmp(:,:) = qatm_ice(:,:)
	411	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
	412	ztmp(:,:)=10.0
	413	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
[2874]	414
	415	! May want to check all values are physically realistic (as in CICE routine
	416	! prepare_forcing)?
	417
	418	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	419	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[2874]	420	CALL nemo2cice(ztmp,swvdr,'T', 1. )
[3193]	421	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[2874]	422	CALL nemo2cice(ztmp,swvdf,'T', 1. )
[3193]	423	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[2874]	424	CALL nemo2cice(ztmp,swidr,'T', 1. )
[3193]	425	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[2874]	426	CALL nemo2cice(ztmp,swidf,'T', 1. )
	427
	428	ENDIF
	429
	430	! Snowfall
[4990]	431	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
	432	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
[3193]	433	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
	434	CALL nemo2cice(ztmp,fsnow,'T', 1. )
[2874]	435
	436	! Rainfall
[4990]	437	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
[3193]	438	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	439	CALL nemo2cice(ztmp,frain,'T', 1. )
[2874]	440
	441	! Freezing/melting potential
[3275]	442	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[6140]	443	nfrzmlt(:,:) = rau0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
[2874]	444
[3193]	445	ztmp(:,:) = nfrzmlt(:,:)
	446	CALL nemo2cice(ztmp,frzmlt,'T', 1. )
[2874]	447
	448	! SST and SSS
	449
[3193]	450	CALL nemo2cice(sst_m,sst,'T', 1. )
	451	CALL nemo2cice(sss_m,sss,'T', 1. )
[2874]	452
	453	! x comp and y comp of surface ocean current
	454	! U point to F point
[3193]	455	DO jj=1,jpjm1
	456	DO ji=1,jpi
	457	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
	458	ENDDO
	459	ENDDO
	460	CALL nemo2cice(ztmp,uocn,'F', -1. )
[2874]	461
	462	! V point to F point
[3193]	463	DO jj=1,jpj
	464	DO ji=1,jpim1
	465	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
	466	ENDDO
	467	ENDDO
	468	CALL nemo2cice(ztmp,vocn,'F', -1. )
[2874]	469
[8313]	470	IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
[3625]	471	!
	472	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	473	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	474	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	475	!
	476	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	477	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	478	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	479	!
	480	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	481	!
	482	!
	483	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	484	zpice(:,:) = ssh_m(:,:)
	485	ENDIF
	486
[3189]	487	! x comp and y comp of sea surface slope (on F points)
	488	! T point to F point
[5836]	489	DO jj = 1, jpjm1
	490	DO ji = 1, jpim1
	491	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) &
	492	& + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1)
	493	END DO
	494	END DO
	495	CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
[3189]	496
	497	! T point to F point
[5836]	498	DO jj = 1, jpjm1
	499	DO ji = 1, jpim1
	500	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) &
	501	& + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1)
	502	END DO
	503	END DO
[3193]	504	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
[3189]	505
[5516]	506	CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
[3152]	507	CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
[3193]	508	!
	509	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_in')
	510	!
[2874]	511	END SUBROUTINE cice_sbc_in
	512
[3152]	513
[5836]	514	SUBROUTINE cice_sbc_out( kt, ksbc )
[2874]	515	!!---------------------------------------------------------------------
	516	!! * ROUTINE cice_sbc_out *
[3040]	517	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
[3152]	518	!!---------------------------------------------------------------------
[2874]	519	INTEGER, INTENT( in ) :: kt ! ocean time step
[4990]	520	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[3152]	521
[3625]	522	INTEGER :: ji, jj, jl ! dummy loop indices
	523	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
[2874]	524	!!---------------------------------------------------------------------
	525
[3193]	526	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_out')
	527	!
[3625]	528	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
[3152]	529
	530	IF( kt == nit000 ) THEN
[2874]	531	IF(lwp) WRITE(numout,*)'cice_sbc_out'
[3152]	532	ENDIF
	533
[2874]	534	! x comp of ocean-ice stress
[3625]	535	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
[3193]	536	ss_iou(:,:)=0.0
[2874]	537	! F point to U point
[3193]	538	DO jj=2,jpjm1
	539	DO ji=2,jpim1
[3625]	540	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
[3193]	541	ENDDO
	542	ENDDO
	543	CALL lbc_lnk( ss_iou , 'U', -1. )
[2874]	544
	545	! y comp of ocean-ice stress
[3625]	546	CALL cice2nemo(strocny,ztmp1,'F', -1. )
[3193]	547	ss_iov(:,:)=0.0
[2874]	548	! F point to V point
	549
[3193]	550	DO jj=1,jpjm1
	551	DO ji=2,jpim1
[3625]	552	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
[3193]	553	ENDDO
	554	ENDDO
	555	CALL lbc_lnk( ss_iov , 'V', -1. )
[2874]	556
	557	! x and y comps of surface stress
	558	! Combine wind stress and ocean-ice stress
	559	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
[5133]	560	! strocnx and strocny already weighted by ice fraction in CICE so not done here
[2874]	561
[3193]	562	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
	563	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
[5133]	564
	565	! Also need ice/ocean stress on T points so that taum can be updated
	566	! This interpolation is already done in CICE so best to use those values
	567	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
	568	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
	569
	570	! Update taum with modulus of ice-ocean stress
	571	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
[5836]	572	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp1ztmp1 + ztmp2ztmp2)
[2874]	573
	574	! Freshwater fluxes
	575
[4990]	576	IF (ksbc == jp_flx) THEN
[2874]	577	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
	578	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
	579	! Not ideal since aice won't be the same as in the atmosphere.
	580	! Better to use evap and tprecip? (but for now don't read in evap in this case)
[3193]	581	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
[7646]	582	ELSE IF (ksbc == jp_blk) THEN
[3193]	583	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
[5407]	584	ELSE IF (ksbc == jp_purecpl) THEN
[3625]	585	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
	586	! This is currently as required with the coupling fields from the UM atmosphere
[3193]	587	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
	588	ENDIF
[2874]	589
[4990]	590	#if defined key_cice4
[3625]	591	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
	592	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
[4990]	593	#else
	594	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
	595	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
	596	#endif
[2874]	597
[3625]	598	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
	599	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
	600	! which has to be compensated for by the ocean salinity potentially going negative
	601	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
	602	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
	603	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
	604	sfx(:,:)=ztmp2(:,:)*1000.0
	605	emp(:,:)=emp(:,:)-ztmp1(:,:)
[4990]	606	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
	607
[3193]	608	CALL lbc_lnk( emp , 'T', 1. )
[3625]	609	CALL lbc_lnk( sfx , 'T', 1. )
[2874]	610
	611	! Solar penetrative radiation and non solar surface heat flux
	612
	613	! Scale qsr and qns according to ice fraction (bulk formulae only)
	614
[7646]	615	IF (ksbc == jp_blk) THEN
[3193]	616	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
	617	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
	618	ENDIF
[2874]	619	! Take into account snow melting except for fully coupled when already in qns_tot
[5407]	620	IF (ksbc == jp_purecpl) THEN
[3193]	621	qsr(:,:)= qsr_tot(:,:)
	622	qns(:,:)= qns_tot(:,:)
	623	ELSE
	624	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
	625	ENDIF
[2874]	626
	627	! Now add in ice / snow related terms
	628	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
[4990]	629	#if defined key_cice4
[3625]	630	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
[4990]	631	#else
	632	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
	633	#endif
[3625]	634	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
[3193]	635	CALL lbc_lnk( qsr , 'T', 1. )
[2874]	636
[3193]	637	DO jj=1,jpj
	638	DO ji=1,jpi
[2874]	639	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
[3193]	640	ENDDO
	641	ENDDO
[2874]	642
[4990]	643	#if defined key_cice4
[3625]	644	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
[4990]	645	#else
	646	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
	647	#endif
[3625]	648	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
[2874]	649
[3193]	650	CALL lbc_lnk( qns , 'T', 1. )
[2874]	651
	652	! Prepare for the following CICE time-step
	653
[3193]	654	CALL cice2nemo(aice,fr_i,'T', 1. )
[5407]	655	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[3625]	656	DO jl=1,ncat
	657	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	658	ENDDO
	659	ENDIF
[2874]	660
	661	! T point to U point
	662	! T point to V point
[3193]	663	DO jj=1,jpjm1
	664	DO ji=1,jpim1
	665	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	666	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	667	ENDDO
	668	ENDDO
[2874]	669
[3193]	670	CALL lbc_lnk ( fr_iu , 'U', 1. )
	671	CALL lbc_lnk ( fr_iv , 'V', 1. )
[2874]	672
[8313]	673	! set the snow+ice mass
	674	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	675	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	676	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	677	snwice_mass_b(:,:) = snwice_mass(:,:)
	678	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
[3625]	679
[2874]	680	! Release work space
	681
[3625]	682	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	683	!
	684	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
	685	!
[2874]	686	END SUBROUTINE cice_sbc_out
	687
[3152]	688
[2874]	689	SUBROUTINE cice_sbc_hadgam( kt )
	690	!!---------------------------------------------------------------------
	691	!! * ROUTINE cice_sbc_hadgam *
[3040]	692	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	693	!!
	694	!!
	695	INTEGER, INTENT( in ) :: kt ! ocean time step
	696	!!---------------------------------------------------------------------
	697
[3625]	698	INTEGER :: jl ! dummy loop index
[3193]	699	INTEGER :: ierror
[2874]	700
[3193]	701	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
	702	!
	703	IF( kt == nit000 ) THEN
[2874]	704	IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
	705	IF( sbc_cpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
[3193]	706	ENDIF
[2874]	707
	708	! ! =========================== !
	709	! ! Prepare Coupling fields !
	710	! ! =========================== !
	711
	712	! x and y comp of ice velocity
	713
[3193]	714	CALL cice2nemo(uvel,u_ice,'F', -1. )
	715	CALL cice2nemo(vvel,v_ice,'F', -1. )
[2874]	716
	717	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	718
	719	! Snow and ice thicknesses (CO_2 and CO_3)
	720
[3625]	721	DO jl = 1,ncat
	722	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
	723	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
[3193]	724	ENDDO
	725	!
	726	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
	727	!
[2874]	728	END SUBROUTINE cice_sbc_hadgam
	729
	730
	731	SUBROUTINE cice_sbc_final
	732	!!---------------------------------------------------------------------
	733	!! * ROUTINE cice_sbc_final *
	734	!! ** Purpose: Finalize CICE
	735	!!---------------------------------------------------------------------
	736
	737	IF(lwp) WRITE(numout,*)'cice_sbc_final'
	738
[3193]	739	CALL CICE_Finalize
[2874]	740
	741	END SUBROUTINE cice_sbc_final
	742
	743	SUBROUTINE cice_sbc_force (kt)
	744	!!---------------------------------------------------------------------
	745	!! * ROUTINE cice_sbc_force *
	746	!! ** Purpose : Provide CICE forcing from files
	747	!!
	748	!!---------------------------------------------------------------------
	749	!! ** Method : READ monthly flux file in NetCDF files
	750	!!
	751	!! snowfall
	752	!! rainfall
	753	!! sublimation rate
	754	!! topmelt (category)
	755	!! botmelt (category)
	756	!!
	757	!! History :
	758	!!----------------------------------------------------------------------
	759	!! * Modules used
	760	USE iom
	761
	762	!! * arguments
	763	INTEGER, INTENT( in ) :: kt ! ocean time step
	764
	765	INTEGER :: ierror ! return error code
	766	INTEGER :: ifpr ! dummy loop index
	767	!!
	768	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	769	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	770	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	771	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	772	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	773
	774	!!
	775	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	776	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	777	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
[4230]	778	INTEGER :: ios
[2874]	779	!!---------------------------------------------------------------------
	780
	781	! ! ====================== !
	782	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	783	! ! ====================== !
[4990]	784	! namsbc_cice is not yet in the reference namelist
	785	! set file information (default values)
	786	cn_dir = './' ! directory in which the model is executed
	787
	788	! (NB: frequency positive => hours, negative => months)
	789	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
	790	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
	791	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	792	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	793	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
	794	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	795	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	796	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	797	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	798	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	799	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	800	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	801	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	802	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	803	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	804
[4230]	805	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
	806	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
	807	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
[2874]	808
[4230]	809	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
	810	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
	811	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
[4624]	812	IF(lwm) WRITE ( numond, namsbc_cice )
[2874]	813
	814	! store namelist information in an array
	815	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	816	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	817	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	818	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	819	slf_i(jp_bot5) = sn_bot5
	820
	821	! set sf structure
	822	ALLOCATE( sf(jpfld), STAT=ierror )
	823	IF( ierror > 0 ) THEN
	824	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	825	ENDIF
	826
	827	DO ifpr= 1, jpfld
	828	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	829	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	830	END DO
	831
	832	! fill sf with slf_i and control print
	833	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	834	!
	835	ENDIF
	836
	837	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	838	! ! input fields at the current time-step
	839
	840	! set the fluxes from read fields
	841	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	842	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	843	! May be better to do this conversion somewhere else
[2874]	844	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	845	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	846	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	847	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	848	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	849	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	850	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	851	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	852	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	853	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	854	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
	855
	856	! control print (if less than 100 time-step asked)
	857	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	858	WRITE(numout,*)
	859	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	860	CALL FLUSH(numout)
	861	ENDIF
	862
	863	END SUBROUTINE cice_sbc_force
	864
	865	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
	866	!!---------------------------------------------------------------------
	867	!! * ROUTINE nemo2cice *
	868	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	869	#if defined key_nemocice_decomp
	870	!!
	871	!! NEMO and CICE PE sub domains are identical, hence
	872	!! there is no need to gather or scatter data from
	873	!! one PE configuration to another.
	874	#else
	875	!! Automatically gather/scatter between
	876	!! different processors and blocks
	877	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	878	!! B. Gather pn into global array (png)
	879	!! C. Map png into CICE global array (pcg)
	880	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	881	#endif
	882	!!---------------------------------------------------------------------
[3193]	883	CHARACTER(len=1), INTENT( in ) :: &
	884	cd_type ! nature of pn grid-point
	885	! ! = T or F gridpoints
	886	REAL(wp), INTENT( in ) :: &
	887	psgn ! control of the sign change
	888	! ! =-1 , the sign is modified following the type of b.c. used
	889	! ! = 1 , no sign change
	890	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	891	#if !defined key_nemocice_decomp
[3625]	892	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	893	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	894	#endif
[3193]	895	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
	896	INTEGER (int_kind) :: &
	897	field_type, &! id for type of field (scalar, vector, angle)
	898	grid_loc ! id for location on horizontal grid
[2874]	899	! (center, NEcorner, Nface, Eface)
	900
[3193]	901	INTEGER :: ji, jj, jn ! dummy loop indices
[5836]	902	!!---------------------------------------------------------------------
[2874]	903
[3193]	904	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	905
[3193]	906	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	907
	908	#if defined key_nemocice_decomp
	909
[3193]	910	! Copy local domain data from NEMO to CICE field
	911	pc(:,:,1)=0.0
[3625]	912	DO jj=2,ny_block-1
	913	DO ji=2,nx_block-1
	914	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	915	ENDDO
	916	ENDDO
[2874]	917
	918	#else
	919
[3193]	920	! B. Gather pn into global array (png)
[2874]	921
[3193]	922	IF ( jpnij > 1) THEN
	923	CALL mppsync
	924	CALL mppgather (pn,0,png)
	925	CALL mppsync
	926	ELSE
	927	png(:,:,1)=pn(:,:)
	928	ENDIF
[2874]	929
[3193]	930	! C. Map png into CICE global array (pcg)
[2874]	931
	932	! Need to make sure this is robust to changes in NEMO halo rows....
	933	! (may be OK but not 100% sure)
	934
[3193]	935	IF (nproc==0) THEN
[2874]	936	! pcg(:,:)=0.0
[3193]	937	DO jn=1,jpnij
[3625]	938	DO jj=nldjt(jn),nlejt(jn)
	939	DO ji=nldit(jn),nleit(jn)
	940	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	941	ENDDO
	942	ENDDO
	943	ENDDO
[3625]	944	DO jj=1,ny_global
	945	DO ji=1,nx_global
	946	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	947	ENDDO
	948	ENDDO
[3193]	949	ENDIF
[2874]	950
	951	#endif
	952
[3193]	953	SELECT CASE ( cd_type )
	954	CASE ( 'T' )
	955	grid_loc=field_loc_center
	956	CASE ( 'F' )
	957	grid_loc=field_loc_NEcorner
	958	END SELECT
[2874]	959
[3193]	960	SELECT CASE ( NINT(psgn) )
	961	CASE ( -1 )
	962	field_type=field_type_vector
	963	CASE ( 1 )
	964	field_type=field_type_scalar
	965	END SELECT
[2874]	966
	967	#if defined key_nemocice_decomp
[3193]	968	! Ensure CICE halos are up to date
	969	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	970	#else
[3193]	971	! D. Scatter pcg to CICE blocks (pc) + update halos
	972	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	973	#endif
	974
	975	END SUBROUTINE nemo2cice
	976
	977	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	978	!!---------------------------------------------------------------------
	979	!! * ROUTINE cice2nemo *
	980	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	981	#if defined key_nemocice_decomp
	982	!!
	983	!! NEMO and CICE PE sub domains are identical, hence
	984	!! there is no need to gather or scatter data from
	985	!! one PE configuration to another.
	986	#else
	987	!! Automatically deal with scatter/gather between
	988	!! different processors and blocks
	989	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	990	!! B. Map pcg into NEMO global array (png)
	991	!! C. Scatter png into NEMO field (pn) for each processor
	992	!! D. Ensure all haloes are filled in pn
	993	#endif
	994	!!---------------------------------------------------------------------
	995
[3193]	996	CHARACTER(len=1), INTENT( in ) :: &
	997	cd_type ! nature of pn grid-point
	998	! ! = T or F gridpoints
	999	REAL(wp), INTENT( in ) :: &
	1000	psgn ! control of the sign change
	1001	! ! =-1 , the sign is modified following the type of b.c. used
	1002	! ! = 1 , no sign change
	1003	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	1004
	1005	#if defined key_nemocice_decomp
[3193]	1006	INTEGER (int_kind) :: &
	1007	field_type, & ! id for type of field (scalar, vector, angle)
	1008	grid_loc ! id for location on horizontal grid
	1009	! (center, NEcorner, Nface, Eface)
[2874]	1010	#else
[3193]	1011	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	1012	#endif
	1013
[3193]	1014	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	1015
[3193]	1016	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	1017
	1018
	1019	#if defined key_nemocice_decomp
	1020
[3193]	1021	SELECT CASE ( cd_type )
	1022	CASE ( 'T' )
	1023	grid_loc=field_loc_center
	1024	CASE ( 'F' )
	1025	grid_loc=field_loc_NEcorner
	1026	END SELECT
[2874]	1027
[3193]	1028	SELECT CASE ( NINT(psgn) )
	1029	CASE ( -1 )
	1030	field_type=field_type_vector
	1031	CASE ( 1 )
	1032	field_type=field_type_scalar
	1033	END SELECT
[2874]	1034
[3193]	1035	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1036
	1037
[3193]	1038	pn(:,:)=0.0
	1039	DO jj=1,jpjm1
	1040	DO ji=1,jpim1
[3625]	1041	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	1042	ENDDO
	1043	ENDDO
[2874]	1044
	1045	#else
	1046
[3193]	1047	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	1048
[3193]	1049	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1050
	1051	! B. Map pcg into NEMO global array (png)
	1052
	1053	! Need to make sure this is robust to changes in NEMO halo rows....
	1054	! (may be OK but not spent much time thinking about it)
[3625]	1055	! Note that non-existent pcg elements may be used below, but
	1056	! the lbclnk call on pn will replace these with sensible values
[2874]	1057
[3193]	1058	IF (nproc==0) THEN
	1059	png(:,:,:)=0.0
	1060	DO jn=1,jpnij
[3625]	1061	DO jj=nldjt(jn),nlejt(jn)
	1062	DO ji=nldit(jn),nleit(jn)
	1063	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1064	ENDDO
	1065	ENDDO
	1066	ENDDO
	1067	ENDIF
[2874]	1068
[3193]	1069	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1070
[3193]	1071	IF ( jpnij > 1) THEN
	1072	CALL mppsync
	1073	CALL mppscatter (png,0,pn)
	1074	CALL mppsync
	1075	ELSE
	1076	pn(:,:)=png(:,:,1)
	1077	ENDIF
[2874]	1078
	1079	#endif
	1080
[3193]	1081	! D. Ensure all haloes are filled in pn
[2874]	1082
[3193]	1083	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1084
	1085	END SUBROUTINE cice2nemo
	1086
	1087	#else
	1088	!!----------------------------------------------------------------------
	1089	!! Default option Dummy module NO CICE sea-ice model
	1090	!!----------------------------------------------------------------------
	1091	CONTAINS
	1092
[4990]	1093	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
[2874]	1094	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1095	END SUBROUTINE sbc_ice_cice
	1096
[4990]	1097	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
[2874]	1098	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1099	END SUBROUTINE cice_sbc_init
	1100
	1101	SUBROUTINE cice_sbc_final ! Dummy routine
	1102	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1103	END SUBROUTINE cice_sbc_final
	1104
	1105	#endif
	1106
	1107	!!======================================================================
	1108	END MODULE sbcice_cice

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

Download in other formats: