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

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

Last change on this file since 3372 was 3372, checked in by sga, 12 years ago

NEMO branch dev_r3337_NOCS10_ICB: change all routine names and create more Gurvanistic havoc

File size: 20.7 KB
Line 
1MODULE 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   !!----------------------------------------------------------------------
45CONTAINS
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 * (njmpp+jj-1) + nimpp+ji-1   !!gm use here mjg(jj) and mig(ji)...
109!!gm        src_calving     (ji,jj) = nicbpack * mjg(jj) + mig(ji)
110         END DO
111      END DO
112      CALL lbc_lnk( src_calving_hflx, 'T', 1._wp )
113      CALL lbc_lnk( src_calving     , 'T', 1._wp )
114
115      ! work out interior of processor from exchange array
116      ! first entry with narea for this processor is left hand interior index
117      ! last  entry                               is right hand interior index
118      jj = jpj/2
119      nicbdi = -1
120      nicbei = -1
121      DO ji = 1, jpi
122         i3 = INT( src_calving(ji,jj) )
123         i2 = INT( i3/nicbpack )
124         i1 = i3 - i2*nicbpack
125         i3 = INT( src_calving_hflx(ji,jj) )
126         IF( i1 == nimpp+ji-1 .AND. i3 == narea ) THEN
127            IF( nicbdi < 0 ) THEN   ;   nicbdi = ji
128            ELSE                    ;   nicbei = ji
129            ENDIF
130         ENDIF
131      END DO
132      !
133      ! repeat for j direction
134      ji = jpi/2
135      nicbdj = -1
136      nicbej = -1
137      DO jj = 1, jpj
138         i3 = INT( src_calving(ji,jj) )
139         i2 = INT( i3/nicbpack )
140         i1 = i3 - i2*nicbpack
141         i3 = INT( src_calving_hflx(ji,jj) )
142         IF( i2 == njmpp+jj-1 .AND. i3 == narea ) THEN
143            IF( nicbdj < 0 ) THEN   ;   nicbdj = jj
144            ELSE                    ;   nicbej = jj
145            ENDIF
146         ENDIF
147      END DO
148      !   
149      ! special for east-west boundary exchange we save the destination index
150      i1 = MAX( nicbdi-1, 1)
151      i3 = INT( src_calving(i1,jpj/2) )
152      jj = INT( i3/nicbpack )
153      ricb_left = REAL( i3 - nicbpack*jj, wp )
154      i1 = MIN( nicbei+1, jpi )
155      i3 = INT( src_calving(i1,jpj/2) )
156      jj = INT( i3/nicbpack )
157      ricb_right = REAL( i3 - nicbpack*jj, wp )
158     
159      ! north fold
160      IF( npolj > 0 ) THEN
161         !
162         ! icebergs in row nicbej+1 get passed across fold
163         nicbfldpts(:)  = INT( src_calving(:,nicbej+1) )
164         nicbflddest(:) = INT( src_calving_hflx(:,nicbej+1) )
165         !
166         ! work out list of unique processors to talk to
167         DO ji = nicbdi, nicbei
168            ii = nicbflddest(ji)
169            DO jn = 1, jpni
170               IF( nicbfldproc(jn) == -1 ) THEN
171                  nicbfldproc(jn) = ii
172                  EXIT                             !!gm EXIT should be avoided: use DO WHILE expression instead
173               ENDIF
174               IF( nicbfldproc(jn) == ii ) EXIT
175            END DO
176         END DO
177      ENDIF
178      !
179      IF( nn_verbose_level > 0) THEN
180         WRITE(numicb,*) 'processor ', narea
181         WRITE(numicb,*) 'jpi, jpj   ', jpi, jpj
182         WRITE(numicb,*) 'nldi, nlei ', nldi, nlei
183         WRITE(numicb,*) 'nldj, nlej ', nldj, nlej
184         WRITE(numicb,*) 'berg i interior ', nicbdi, nicbei
185         WRITE(numicb,*) 'berg j interior ', nicbdj, nicbej
186         WRITE(numicb,*) 'berg left       ', ricb_left
187         WRITE(numicb,*) 'berg right      ', ricb_right
188         jj = jpj/2
189         WRITE(numicb,*) "central j line:"
190         WRITE(numicb,*) "i processor"
191         WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), ji=1,jpi)
192         WRITE(numicb,*) "i point"
193         WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi)
194         ji = jpi/2
195         WRITE(numicb,*) "central i line:"
196         WRITE(numicb,*) "j processor"
197         WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), jj=1,jpj)
198         WRITE(numicb,*) "j point"
199         WRITE(numicb,*) (INT(src_calving(ji,jj)), jj=1,jpj)
200         IF( npolj > 0 ) THEN
201            WRITE(numicb,*) 'north fold destination points '
202            WRITE(numicb,*) nicbfldpts
203            WRITE(numicb,*) 'north fold destination procs  '
204            WRITE(numicb,*) nicbflddest
205         ENDIF
206         CALL flush(numicb)
207      ENDIF
208     
209      src_calving     (:,:) = 0._wp
210      src_calving_hflx(:,:) = 0._wp
211
212      ! assign each new iceberg with a unique number constructed from the processor number
213      ! and incremented by the total number of processors
214      num_bergs(:) = 0
215      num_bergs(1) = narea - jpnij
216
217      ! when not generating test icebergs we need to setup calving file
218      IF( nn_test_icebergs < 0 ) THEN
219         !
220         ! maximum distribution class array does not change in time so read it once
221         cl_sdist = TRIM( cn_dir )//TRIM( sn_icb%clname )
222         CALL iom_open ( cl_sdist, inum )                              ! open file
223         ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. )
224         IF( ivar > 0 ) THEN
225            CALL iom_get  ( inum, jpdom_data, 'maxclass', src_calving )   ! read the max distribution array
226            berg_grid%maxclass(:,:) = INT( src_calving )
227            src_calving(:,:) = 0._wp
228         ENDIF
229         CALL iom_close( inum )                                     ! close file
230         !
231         WRITE(numicb,*)
232         WRITE(numicb,*) '          calving read in a file'
233         ALLOCATE( sf_icb(1), STAT=istat1 )         ! Create sf_icb structure (calving)
234         ALLOCATE( sf_icb(1)%fnow(jpi,jpj,1), STAT=istat2 )
235         ALLOCATE( sf_icb(1)%fdta(jpi,jpj,1,2), STAT=istat3 )
236         IF( istat1+istat2+istat3 > 0 ) THEN
237            CALL ctl_stop( 'sbc_icb: unable to allocate sf_icb structure' )   ;   RETURN
238         ENDIF
239         !                                          ! fill sf_icb with the namelist (sn_icb) and control print
240         CALL fld_fill( sf_icb, (/ sn_icb /), cn_dir, 'icb_init', 'read calving data', 'namicb' )
241         !
242      ENDIF
243
244      IF( .NOT.ln_rstart ) THEN
245         IF( nn_test_icebergs > 0 )   CALL icb_ini_gen()
246      ELSE
247         IF( nn_test_icebergs > 0 ) THEN
248            CALL icb_ini_gen()
249         ELSE
250            CALL icb_rst_read()
251            l_restarted_bergs = .TRUE.
252         ENDIF
253      ENDIF
254      !
255      IF( nn_sample_rate .GT. 0 ) CALL icb_trj_init( nitend )
256      !
257      CALL icb_dia_init()
258      !
259      IF( nn_verbose_level >= 2 )   CALL icb_utl_print('icb_init, initial status', nit000-1)
260      !
261   END SUBROUTINE icb_init
262
263   SUBROUTINE icb_ini_gen()
264      !!----------------------------------------------------------------------
265      !!                  ***  ROUTINE icb_ini_gen  ***
266      !!
267      !! ** Purpose :   iceberg generation
268      !!
269      !! ** Method  : - at each grid point of the test box supplied in the namelist
270      !!                generate an iceberg in one class determined by the value of
271      !!                parameter nn_test_icebergs
272      !!----------------------------------------------------------------------
273      INTEGER                         ::   ji, jj, ibergs
274      TYPE(iceberg)                   ::   localberg ! NOT a pointer but an actual local variable
275      TYPE(point)                     ::   localpt
276      INTEGER                         ::   iyr, imon, iday, ihr, imin, isec
277      INTEGER                         ::   iberg
278      !!----------------------------------------------------------------------
279
280      ! For convenience
281      iberg = nn_test_icebergs
282
283      ! call get_date(Time, iyr, imon, iday, ihr, imin, isec)
284      ! Convert nemo time variables from dom_oce into local versions
285      iyr  = nyear
286      imon = nmonth
287      iday = nday
288      ihr = INT(nsec_day/3600)
289      imin = INT((nsec_day-ihr*3600)/60)
290      isec = nsec_day - ihr*3600 - imin*60
291
292      ! no overlap for icebergs since we want only one instance of each across the whole domain
293      ! so restrict area of interest
294      ! use tmask here because tmask_i has been doctored on one side of the north fold line
295
296      DO jj = nicbdj, nicbej
297         DO ji = nicbdi, nicbei
298            IF( tmask(ji,jj,1) > 0._wp        .AND.                                       &
299                rn_test_box(1) < glamt(ji,jj) .AND. glamt(ji,jj) < rn_test_box(2) .AND.   &
300                rn_test_box(3) < gphit(ji,jj) .AND. gphit(ji,jj) < rn_test_box(4) ) THEN
301               localberg%mass_scaling = rn_mass_scaling(iberg)
302               localpt%xi = REAL( nimpp+ji-1, wp )
303               localpt%yj = REAL( njmpp+jj-1, wp )
304               localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T', 0, 0 )
305               localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T', 0, 0 )
306               localpt%mass      = rn_initial_mass     (iberg)
307               localpt%thickness = rn_initial_thickness(iberg)
308               localpt%width  = first_width (iberg)
309               localpt%length = first_length(iberg)
310               localpt%year = iyr
311               localpt%day = FLOAT(iday)+(FLOAT(ihr)+FLOAT(imin)/60._wp)/24._wp
312               localpt%mass_of_bits = 0._wp
313               localpt%heat_density = 0._wp
314               localpt%uvel = 0._wp
315               localpt%vvel = 0._wp
316               CALL icb_utl_incr()
317               localberg%number(:) = num_bergs(:)
318               call icb_utl_add(localberg, localpt)
319            ENDIF
320         END DO
321      END DO
322      !
323      ibergs = icb_utl_count()
324      IF( lk_mpp ) CALL mpp_sum(ibergs)
325      WRITE(numicb,'(a,i6,a)') 'diamonds, icb_ini_gen: ',ibergs,' were generated'
326      !
327   END SUBROUTINE icb_ini_gen
328
329   SUBROUTINE icb_nam
330      !!----------------------------------------------------------------------
331      !!                     ***  ROUTINE icb_nam  ***
332      !!
333      !! ** Purpose :   read iceberg namelist and print the variables.
334      !!
335      !! ** input   : - namberg namelist
336      !!----------------------------------------------------------------------
337      INTEGER  ::   jn      ! dummy loop indices
338      REAL(wp) ::   zfact   ! local scalar
339      !
340      NAMELIST/namberg/ ln_icebergs    , ln_bergdia     , nn_sample_rate      , rn_initial_mass      ,   &
341         &              rn_distribution, rn_mass_scaling, rn_initial_thickness, nn_verbose_write     ,   &
342         &              rn_rho_bergs   , rn_LoW_ratio   , nn_verbose_level    , ln_operator_splitting,   &
343         &              rn_bits_erosion_fraction        , rn_sicn_shift       , ln_passive_mode      ,   &
344         &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   &
345         &              rn_speed_limit , cn_dir, sn_icb
346      !!----------------------------------------------------------------------
347
348      ! (NB: frequency positive => hours, negative => months)
349      !            !   file     ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   !
350      !            !   name     !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      !
351      sn_icb = FLD_N( 'calving' ,    -1     , 'calving'  ,  .TRUE.    , .TRUE. ,   'yearly'  , ''       , ''         )
352
353      REWIND( numnam )              ! Namelist namrun : iceberg parameters
354      READ  ( numnam, namberg )
355     
356!!gm  bug here in mpp : the return instruction will be done only on the 1st processor (lwp=True)
357!      IF( .NOT. ln_icebergs .AND. lwp ) THEN
358!         WRITE(numout,*)
359!         WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
360!         WRITE(numout,*) '~~~~~~~~ '
361!         WRITE(numout,*) 'NO icebergs used'
362!         RETURN
363!      ENDIF
364!!gm solution:
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!!gm end
375
376      IF( nn_test_icebergs > nclasses ) THEN
377          IF(lwp) WRITE(numout,*) 'Resetting nn_test_icebergs to ', nclasses
378          nn_test_icebergs = nclasses
379      ENDIF
380
381      zfact = SUM( rn_distribution )
382      IF( zfact < 1._wp ) THEN
383         IF( zfact <= 0._wp ) THEN
384            CALL ctl_stop( 'icb_nam: sum of berg distribution equal to zero' )
385         ELSE
386            rn_distribution(:) = rn_distribution(:) / zfact
387            CALL ctl_warn( 'icb_nam: sum of berg input distribution not equal to one and so RESCALED' )
388         ENDIF
389      ENDIF
390
391      IF(lwp) THEN                  ! control print
392         WRITE(numout,*)
393         WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
394         WRITE(numout,*) '~~~~~~~~ '
395         WRITE(numout,*) '   Calculate budgets                                            ln_bergdia       = ', ln_bergdia
396         WRITE(numout,*) '   Period between sampling of position for trajectory storage   nn_sample_rate = ', nn_sample_rate
397         WRITE(numout,*) '   Mass thresholds between iceberg classes (kg)                 rn_initial_mass     ='
398         DO jn=1,nclasses
399            WRITE(numout,'(a,f15.2)') '                                                                ',rn_initial_mass(jn)
400         ENDDO
401         WRITE(numout,*) '   Fraction of calving to apply to this class (non-dim)         rn_distribution     ='
402         DO jn = 1, nclasses
403            WRITE(numout,'(a,f10.2)') '                                                                ',rn_distribution(jn)
404         END DO
405         WRITE(numout,*) '   Ratio between effective and real iceberg mass (non-dim)      rn_mass_scaling     = '
406         DO jn = 1, nclasses
407            WRITE(numout,'(a,f10.2)') '                                                                ',rn_mass_scaling(jn)
408         END DO
409         WRITE(numout,*) '   Total thickness of newly calved bergs (m)                    rn_initial_thickness = '
410         DO jn = 1, nclasses
411            WRITE(numout,'(a,f10.2)') '                                                                ',rn_initial_thickness(jn)
412         END DO
413         WRITE(numout,*) '   Timesteps between verbose messages                           nn_verbose_write    = ', nn_verbose_write
414
415         WRITE(numout,*) '   Density of icebergs                           rn_rho_bergs  = ', rn_rho_bergs
416         WRITE(numout,*) '   Initial ratio L/W for newly calved icebergs   rn_LoW_ratio  = ', rn_LoW_ratio
417         WRITE(numout,*) '   Turn on more verbose output                          level  = ', nn_verbose_level
418         WRITE(numout,*) '   Use first order operator splitting for thermodynamics    ',   &
419            &                    'use_operator_splitting = ', ln_operator_splitting
420         WRITE(numout,*) '   Fraction of erosion melt flux to divert to bergy bits    ',   &
421            &                    'bits_erosion_fraction = ', rn_bits_erosion_fraction
422
423         WRITE(numout,*) '   Shift of sea-ice concentration in erosion flux modulation ',   &
424            &                    '(0<sicn_shift<1)    rn_sicn_shift  = ', rn_sicn_shift
425         WRITE(numout,*) '   Do not add freshwater flux from icebergs to ocean                ',   &
426            &                    '                  passive_mode            = ', ln_passive_mode
427         WRITE(numout,*) '   Time average the weight on the ocean   time_average_weight       = ', ln_time_average_weight
428         WRITE(numout,*) '   Create icebergs in absence of a restart file   nn_test_icebergs  = ', nn_test_icebergs
429         WRITE(numout,*) '                   in lon/lat box                                   = ', rn_test_box
430         WRITE(numout,*) '   CFL speed limit for a berg            speed_limit                = ', rn_speed_limit
431         WRITE(numout,*) '   Writing Iceberg status information to icebergs.stat file        '
432      ENDIF
433      !
434   END SUBROUTINE icb_nam
435
436   !!======================================================================
437
438END MODULE icbini
Note: See TracBrowser for help on using the repository browser.