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.

icbini.F90 in branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

Context Navigation

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90 @ 3374

Last change on this file since 3374 was 3374, checked in by sga, 12 years ago
NEMO branch dev_r3337_NOCS10_ICB: adopt mig, mjg, mi1 and mj1 for improved obscurity
File size: 20.6 KB

Line
1	MODULE icbini
2
3	!!======================================================================
4	!! * MODULE icbini *
5	!! Icebergs: initialise variables for iceberg tracking
6	!!======================================================================
7	!! History : - ! 2010-01 (T. Martin & A. Adcroft) Original code
8	!! 3.3 ! 2011-03 (G. Madec) Part conversion to NEMO form ; Removal of mapping from another grid
9	!! - ! 2011-04 (S. Alderson) Split into separate modules ; Restore restart routines
10	!! - ! 2011-05 (S. Alderson) generate_test_icebergs restored ; new forcing arrays with extra halo ;
11	!! - ! north fold exchange arrays added
12	!!----------------------------------------------------------------------
13	!!----------------------------------------------------------------------
14	!! icb_init : initialise icebergs
15	!! icb_ini_gen : generate test icebergs
16	!! icb_nam : read iceberg namelist
17	!!----------------------------------------------------------------------
18	USE dom_oce ! ocean domain
19	USE in_out_manager ! IO routines and numout in particular
20	USE lib_mpp ! mpi library and lk_mpp in particular
21	USE sbc_oce ! ocean: surface boundary condition
22	USE iom ! IOM library
23	USE fldread ! field read
24	USE lbclnk ! lateral boundary condition - MPP link
25	!
26	USE icb_oce ! define iceberg arrays
27	USE icbutl ! iceberg utility routines
28	USE icbrst ! iceberg restart routines
29	USE icbtrj ! iceberg trajectory I/O routines
30	USE icbdia ! iceberg budget routines
31
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC icb_init ! routine called in nemogcm.F90 module
36
37	CHARACTER(len=100) :: cn_dir = './' ! Root directory for location of icb files
38	TYPE(FLD_N) :: sn_icb ! information about the calving file to be read
39
40	!!----------------------------------------------------------------------
41	!! NEMO/OPA 3.3 , NEMO Consortium (2011)
42	!! $Id:$
43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44	!!----------------------------------------------------------------------
45	CONTAINS
46
47	SUBROUTINE icb_init( pdt, kt )
48	!!----------------------------------------------------------------------
49	!! * ROUTINE dom_init *
50	!!
51	!! ** Purpose : iceberg initialization.
52	!!
53	!! ** Method : - blah blah
54	!!----------------------------------------------------------------------
55	REAL(wp), INTENT(in) :: pdt ! iceberg time-step (rdt*nn_fsbc)
56	INTEGER , INTENT(in) :: kt ! time step number
57	!
58	INTEGER :: ji, jj, jn ! dummy loop indices
59	INTEGER :: i1, i2, i3 ! local integers
60	INTEGER :: ii, inum, ivar ! - -
61	INTEGER :: istat1, istat2, istat3 ! - -
62	CHARACTER(len=300) :: cl_sdist ! local character
63	!!----------------------------------------------------------------------
64	!
65	CALL icb_nam ! Read and print namelist parameters
66	!
67	IF( .NOT. ln_icebergs ) RETURN
68
69	! ! allocate gridded fields
70	IF( icb_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' )
71
72	! ! open ascii output file or files for iceberg status information
73	! ! note that we choose to do this on all processors since we cannot
74	! ! predict where icebergs will be ahead of time
75	CALL ctl_opn( numicb, 'icebergs.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
76
77	! set parameters (mostly from namelist)
78	!
79	berg_dt = pdt
80	first_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) ) )
81	first_length(:) = rn_LoW_ratio * first_width(:)
82
83	berg_grid%calving (:,:) = 0._wp
84	berg_grid%calving_hflx (:,:) = 0._wp
85	berg_grid%stored_heat (:,:) = 0._wp
86	berg_grid%floating_melt(:,:) = 0._wp
87	berg_grid%maxclass (:,:) = nclasses
88	berg_grid%stored_ice (:,:,:) = 0._wp
89	berg_grid%tmp (:,:) = 0._wp
90	src_calving (:,:) = 0._wp
91	src_calving_hflx (:,:) = 0._wp
92
93	! ! domain for icebergs
94	IF( lk_mpp .AND. jpni == 1 ) CALL ctl_stop( 'icbinit: having ONE processor in x currently does not work' )
95	! NB: the issue here is simply that cyclic east-west boundary condition have not been coded in mpp case
96	! for the north fold we work out which points communicate by asking
97	! lbc_lnk to pass processor number (valid even in single processor case)
98	! borrow src_calving arrays for this
99	!
100	! pack i and j together using a scaling of a power of 10
101	nicbpack = 10000
102	IF( jpiglo >= nicbpack ) CALL ctl_stop( 'icbini: processor index packing failure' )
103	nicbfldproc(:) = -1
104
105	DO jj = 1, jpj
106	DO ji = 1, jpi
107	src_calving_hflx(ji,jj) = narea
108	src_calving (ji,jj) = nicbpack * mjg(jj) + mig(ji)
109	END DO
110	END DO
111	CALL lbc_lnk( src_calving_hflx, 'T', 1._wp )
112	CALL lbc_lnk( src_calving , 'T', 1._wp )
113
114	! work out interior of processor from exchange array
115	! first entry with narea for this processor is left hand interior index
116	! last entry is right hand interior index
117	jj = jpj/2
118	nicbdi = -1
119	nicbei = -1
120	DO ji = 1, jpi
121	i3 = INT( src_calving(ji,jj) )
122	i2 = INT( i3/nicbpack )
123	i1 = i3 - i2*nicbpack
124	i3 = INT( src_calving_hflx(ji,jj) )
125	IF( i1 == mig(ji) .AND. i3 == narea ) THEN
126	IF( nicbdi < 0 ) THEN ; nicbdi = ji
127	ELSE ; nicbei = ji
128	ENDIF
129	ENDIF
130	END DO
131	!
132	! repeat for j direction
133	ji = jpi/2
134	nicbdj = -1
135	nicbej = -1
136	DO jj = 1, jpj
137	i3 = INT( src_calving(ji,jj) )
138	i2 = INT( i3/nicbpack )
139	i1 = i3 - i2*nicbpack
140	i3 = INT( src_calving_hflx(ji,jj) )
141	IF( i2 == mjg(jj) .AND. i3 == narea ) THEN
142	IF( nicbdj < 0 ) THEN ; nicbdj = jj
143	ELSE ; nicbej = jj
144	ENDIF
145	ENDIF
146	END DO
147	!
148	! special for east-west boundary exchange we save the destination index
149	i1 = MAX( nicbdi-1, 1)
150	i3 = INT( src_calving(i1,jpj/2) )
151	jj = INT( i3/nicbpack )
152	ricb_left = REAL( i3 - nicbpack*jj, wp )
153	i1 = MIN( nicbei+1, jpi )
154	i3 = INT( src_calving(i1,jpj/2) )
155	jj = INT( i3/nicbpack )
156	ricb_right = REAL( i3 - nicbpack*jj, wp )
157
158	! north fold
159	IF( npolj > 0 ) THEN
160	!
161	! icebergs in row nicbej+1 get passed across fold
162	nicbfldpts(:) = INT( src_calving(:,nicbej+1) )
163	nicbflddest(:) = INT( src_calving_hflx(:,nicbej+1) )
164	!
165	! work out list of unique processors to talk to
166	DO ji = nicbdi, nicbei
167	ii = nicbflddest(ji)
168	DO jn = 1, jpni
169	IF( nicbfldproc(jn) == -1 ) THEN
170	nicbfldproc(jn) = ii
171	EXIT !!gm EXIT should be avoided: use DO WHILE expression instead
172	ENDIF
173	IF( nicbfldproc(jn) == ii ) EXIT
174	END DO
175	END DO
176	ENDIF
177	!
178	IF( nn_verbose_level > 0) THEN
179	WRITE(numicb,*) 'processor ', narea
180	WRITE(numicb,*) 'jpi, jpj ', jpi, jpj
181	WRITE(numicb,*) 'nldi, nlei ', nldi, nlei
182	WRITE(numicb,*) 'nldj, nlej ', nldj, nlej
183	WRITE(numicb,*) 'berg i interior ', nicbdi, nicbei
184	WRITE(numicb,*) 'berg j interior ', nicbdj, nicbej
185	WRITE(numicb,*) 'berg left ', ricb_left
186	WRITE(numicb,*) 'berg right ', ricb_right
187	jj = jpj/2
188	WRITE(numicb,*) "central j line:"
189	WRITE(numicb,*) "i processor"
190	WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), ji=1,jpi)
191	WRITE(numicb,*) "i point"
192	WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi)
193	ji = jpi/2
194	WRITE(numicb,*) "central i line:"
195	WRITE(numicb,*) "j processor"
196	WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), jj=1,jpj)
197	WRITE(numicb,*) "j point"
198	WRITE(numicb,*) (INT(src_calving(ji,jj)), jj=1,jpj)
199	IF( npolj > 0 ) THEN
200	WRITE(numicb,*) 'north fold destination points '
201	WRITE(numicb,*) nicbfldpts
202	WRITE(numicb,*) 'north fold destination procs '
203	WRITE(numicb,*) nicbflddest
204	ENDIF
205	CALL flush(numicb)
206	ENDIF
207
208	src_calving (:,:) = 0._wp
209	src_calving_hflx(:,:) = 0._wp
210
211	! assign each new iceberg with a unique number constructed from the processor number
212	! and incremented by the total number of processors
213	num_bergs(:) = 0
214	num_bergs(1) = narea - jpnij
215
216	! when not generating test icebergs we need to setup calving file
217	IF( nn_test_icebergs < 0 ) THEN
218	!
219	! maximum distribution class array does not change in time so read it once
220	cl_sdist = TRIM( cn_dir )//TRIM( sn_icb%clname )
221	CALL iom_open ( cl_sdist, inum ) ! open file
222	ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. )
223	IF( ivar > 0 ) THEN
224	CALL iom_get ( inum, jpdom_data, 'maxclass', src_calving ) ! read the max distribution array
225	berg_grid%maxclass(:,:) = INT( src_calving )
226	src_calving(:,:) = 0._wp
227	ENDIF
228	CALL iom_close( inum ) ! close file
229	!
230	WRITE(numicb,*)
231	WRITE(numicb,*) ' calving read in a file'
232	ALLOCATE( sf_icb(1), STAT=istat1 ) ! Create sf_icb structure (calving)
233	ALLOCATE( sf_icb(1)%fnow(jpi,jpj,1), STAT=istat2 )
234	ALLOCATE( sf_icb(1)%fdta(jpi,jpj,1,2), STAT=istat3 )
235	IF( istat1+istat2+istat3 > 0 ) THEN
236	CALL ctl_stop( 'sbc_icb: unable to allocate sf_icb structure' ) ; RETURN
237	ENDIF
238	! ! fill sf_icb with the namelist (sn_icb) and control print
239	CALL fld_fill( sf_icb, (/ sn_icb /), cn_dir, 'icb_init', 'read calving data', 'namicb' )
240	!
241	ENDIF
242
243	IF( .NOT.ln_rstart ) THEN
244	IF( nn_test_icebergs > 0 ) CALL icb_ini_gen()
245	ELSE
246	IF( nn_test_icebergs > 0 ) THEN
247	CALL icb_ini_gen()
248	ELSE
249	CALL icb_rst_read()
250	l_restarted_bergs = .TRUE.
251	ENDIF
252	ENDIF
253	!
254	IF( nn_sample_rate .GT. 0 ) CALL icb_trj_init( nitend )
255	!
256	CALL icb_dia_init()
257	!
258	IF( nn_verbose_level >= 2 ) CALL icb_utl_print('icb_init, initial status', nit000-1)
259	!
260	END SUBROUTINE icb_init
261
262	SUBROUTINE icb_ini_gen()
263	!!----------------------------------------------------------------------
264	!! * ROUTINE icb_ini_gen *
265	!!
266	!! ** Purpose : iceberg generation
267	!!
268	!! ** Method : - at each grid point of the test box supplied in the namelist
269	!! generate an iceberg in one class determined by the value of
270	!! parameter nn_test_icebergs
271	!!----------------------------------------------------------------------
272	INTEGER :: ji, jj, ibergs
273	TYPE(iceberg) :: localberg ! NOT a pointer but an actual local variable
274	TYPE(point) :: localpt
275	INTEGER :: iyr, imon, iday, ihr, imin, isec
276	INTEGER :: iberg
277	!!----------------------------------------------------------------------
278
279	! For convenience
280	iberg = nn_test_icebergs
281
282	! call get_date(Time, iyr, imon, iday, ihr, imin, isec)
283	! Convert nemo time variables from dom_oce into local versions
284	iyr = nyear
285	imon = nmonth
286	iday = nday
287	ihr = INT(nsec_day/3600)
288	imin = INT((nsec_day-ihr*3600)/60)
289	isec = nsec_day - ihr3600 - imin60
290
291	! no overlap for icebergs since we want only one instance of each across the whole domain
292	! so restrict area of interest
293	! use tmask here because tmask_i has been doctored on one side of the north fold line
294
295	DO jj = nicbdj, nicbej
296	DO ji = nicbdi, nicbei
297	IF( tmask(ji,jj,1) > 0._wp .AND. &
298	rn_test_box(1) < glamt(ji,jj) .AND. glamt(ji,jj) < rn_test_box(2) .AND. &
299	rn_test_box(3) < gphit(ji,jj) .AND. gphit(ji,jj) < rn_test_box(4) ) THEN
300	localberg%mass_scaling = rn_mass_scaling(iberg)
301	localpt%xi = REAL( mig(ji), wp )
302	localpt%yj = REAL( mjg(jj), wp )
303	localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T', 0, 0 )
304	localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T', 0, 0 )
305	localpt%mass = rn_initial_mass (iberg)
306	localpt%thickness = rn_initial_thickness(iberg)
307	localpt%width = first_width (iberg)
308	localpt%length = first_length(iberg)
309	localpt%year = iyr
310	localpt%day = FLOAT(iday)+(FLOAT(ihr)+FLOAT(imin)/60._wp)/24._wp
311	localpt%mass_of_bits = 0._wp
312	localpt%heat_density = 0._wp
313	localpt%uvel = 0._wp
314	localpt%vvel = 0._wp
315	CALL icb_utl_incr()
316	localberg%number(:) = num_bergs(:)
317	call icb_utl_add(localberg, localpt)
318	ENDIF
319	END DO
320	END DO
321	!
322	ibergs = icb_utl_count()
323	IF( lk_mpp ) CALL mpp_sum(ibergs)
324	WRITE(numicb,'(a,i6,a)') 'diamonds, icb_ini_gen: ',ibergs,' were generated'
325	!
326	END SUBROUTINE icb_ini_gen
327
328	SUBROUTINE icb_nam
329	!!----------------------------------------------------------------------
330	!! * ROUTINE icb_nam *
331	!!
332	!! ** Purpose : read iceberg namelist and print the variables.
333	!!
334	!! ** input : - namberg namelist
335	!!----------------------------------------------------------------------
336	INTEGER :: jn ! dummy loop indices
337	REAL(wp) :: zfact ! local scalar
338	!
339	NAMELIST/namberg/ ln_icebergs , ln_bergdia , nn_sample_rate , rn_initial_mass , &
340	& rn_distribution, rn_mass_scaling, rn_initial_thickness, nn_verbose_write , &
341	& rn_rho_bergs , rn_LoW_ratio , nn_verbose_level , ln_operator_splitting, &
342	& rn_bits_erosion_fraction , rn_sicn_shift , ln_passive_mode , &
343	& ln_time_average_weight , nn_test_icebergs , rn_test_box , &
344	& rn_speed_limit , cn_dir, sn_icb
345	!!----------------------------------------------------------------------
346
347	! (NB: frequency positive => hours, negative => months)
348	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation !
349	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs !
350	sn_icb = FLD_N( 'calving' , -1 , 'calving' , .TRUE. , .TRUE. , 'yearly' , '' , '' )
351
352	REWIND( numnam ) ! Namelist namrun : iceberg parameters
353	READ ( numnam, namberg )
354
355	!!gm bug here in mpp : the return instruction will be done only on the 1st processor (lwp=True)
356	! IF( .NOT. ln_icebergs .AND. lwp ) THEN
357	! WRITE(numout,*)
358	! WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
359	! WRITE(numout,*) '~~~~~~~~ '
360	! WRITE(numout,*) 'NO icebergs used'
361	! RETURN
362	! ENDIF
363	!!gm solution:
364
365	IF( .NOT. ln_icebergs ) THEN ! no icebergs
366	IF(lwp) THEN
367	WRITE(numout,*)
368	WRITE(numout,*) 'icb_nam : ln_icebergs = F , NO icebergs used'
369	WRITE(numout,*) '~~~~~~~~ '
370	ENDIF
371	RETURN
372	ENDIF
373	!!gm end
374
375	IF( nn_test_icebergs > nclasses ) THEN
376	IF(lwp) WRITE(numout,*) 'Resetting nn_test_icebergs to ', nclasses
377	nn_test_icebergs = nclasses
378	ENDIF
379
380	zfact = SUM( rn_distribution )
381	IF( zfact < 1._wp ) THEN
382	IF( zfact <= 0._wp ) THEN
383	CALL ctl_stop( 'icb_nam: sum of berg distribution equal to zero' )
384	ELSE
385	rn_distribution(:) = rn_distribution(:) / zfact
386	CALL ctl_warn( 'icb_nam: sum of berg input distribution not equal to one and so RESCALED' )
387	ENDIF
388	ENDIF
389
390	IF(lwp) THEN ! control print
391	WRITE(numout,*)
392	WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
393	WRITE(numout,*) '~~~~~~~~ '
394	WRITE(numout,*) ' Calculate budgets ln_bergdia = ', ln_bergdia
395	WRITE(numout,*) ' Period between sampling of position for trajectory storage nn_sample_rate = ', nn_sample_rate
396	WRITE(numout,*) ' Mass thresholds between iceberg classes (kg) rn_initial_mass ='
397	DO jn=1,nclasses
398	WRITE(numout,'(a,f15.2)') ' ',rn_initial_mass(jn)
399	ENDDO
400	WRITE(numout,*) ' Fraction of calving to apply to this class (non-dim) rn_distribution ='
401	DO jn = 1, nclasses
402	WRITE(numout,'(a,f10.2)') ' ',rn_distribution(jn)
403	END DO
404	WRITE(numout,*) ' Ratio between effective and real iceberg mass (non-dim) rn_mass_scaling = '
405	DO jn = 1, nclasses
406	WRITE(numout,'(a,f10.2)') ' ',rn_mass_scaling(jn)
407	END DO
408	WRITE(numout,*) ' Total thickness of newly calved bergs (m) rn_initial_thickness = '
409	DO jn = 1, nclasses
410	WRITE(numout,'(a,f10.2)') ' ',rn_initial_thickness(jn)
411	END DO
412	WRITE(numout,*) ' Timesteps between verbose messages nn_verbose_write = ', nn_verbose_write
413
414	WRITE(numout,*) ' Density of icebergs rn_rho_bergs = ', rn_rho_bergs
415	WRITE(numout,*) ' Initial ratio L/W for newly calved icebergs rn_LoW_ratio = ', rn_LoW_ratio
416	WRITE(numout,*) ' Turn on more verbose output level = ', nn_verbose_level
417	WRITE(numout,*) ' Use first order operator splitting for thermodynamics ', &
418	& 'use_operator_splitting = ', ln_operator_splitting
419	WRITE(numout,*) ' Fraction of erosion melt flux to divert to bergy bits ', &
420	& 'bits_erosion_fraction = ', rn_bits_erosion_fraction
421
422	WRITE(numout,*) ' Shift of sea-ice concentration in erosion flux modulation ', &
423	& '(0<sicn_shift<1) rn_sicn_shift = ', rn_sicn_shift
424	WRITE(numout,*) ' Do not add freshwater flux from icebergs to ocean ', &
425	& ' passive_mode = ', ln_passive_mode
426	WRITE(numout,*) ' Time average the weight on the ocean time_average_weight = ', ln_time_average_weight
427	WRITE(numout,*) ' Create icebergs in absence of a restart file nn_test_icebergs = ', nn_test_icebergs
428	WRITE(numout,*) ' in lon/lat box = ', rn_test_box
429	WRITE(numout,*) ' CFL speed limit for a berg speed_limit = ', rn_speed_limit
430	WRITE(numout,*) ' Writing Iceberg status information to icebergs.stat file '
431	ENDIF
432	!
433	END SUBROUTINE icb_nam
434
435	!!======================================================================
436
437	END MODULE icbini

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

Download in other formats: