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/UKMO/dev_r5518_CICE_coupling_GSI7_GSI8/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_CICE_coupling_GSI7_GSI8/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5677

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

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

Download in other formats: