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/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

source: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90 @ 4733

Last change on this file since 4733 was 4624, checked in by acc, 10 years ago

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