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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90 @ 3983

Last change on this file since 3983 was 3983, checked in by acc, 11 years ago

Update xml files with missing iceberg (ICB) variables and add new axis definition to iom.F90. The latter requires changes in icb modules to avoid a cyclic dependency

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