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

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago
#2050 first version. Compiled OK in moci test suite
File size: 48.6 KB

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

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

Download in other formats: