MODULE icbini !!====================================================================== !! *** MODULE icbini *** !! Ocean physics: initialise variables for iceberg tracking !!====================================================================== !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code !! - ! 2011-03 (Madec) Part conversion to NEMO form !! - ! Removal of mapping from another grid !! - ! 2011-04 (Alderson) Split into separate modules !! - ! Restore restart routines !! - ! 2011-05 (Alderson) generate_test_icebergs restored !! - ! 2011-05 (Alderson) new forcing arrays with extra halo !! - ! 2011-05 (Alderson) north fold exchange arrays added !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! icb_init : initialise !! icb_gen : generate test icebergs !! icb_nam : read iceberg namelist !!---------------------------------------------------------------------- USE dom_oce ! ocean domain USE in_out_manager ! IO routines and numout in particular USE lib_mpp ! mpi library and lk_mpp in particular USE sbc_oce USE iom USE fldread USE lbclnk USE icb_oce ! define iceberg arrays USE icbutl ! iceberg utility routines USE icbrst ! iceberg restart routines USE icbtrj ! iceberg trajectory I/O routines USE icbdia ! iceberg budget routines IMPLICIT NONE PRIVATE CHARACTER(len=100), PRIVATE :: cn_dir = './' ! Root directory for location of icb files TYPE(FLD_N), PRIVATE :: sn_icb ! information about the calving file to be read PUBLIC icb_init ! routine called in nemogcm.F90 module PUBLIC icb_gen ! routine called in icbclv.F90 module PRIVATE icb_nam PRIVATE icb_alloc CONTAINS INTEGER FUNCTION icb_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE icb_alloc *** !!---------------------------------------------------------------------- ! INTEGER :: ill ! ! open ascii output file or files for iceberg status information ! note that we choose to do this on all processors since we cannot predict where ! icebergs will be ahead of time ! CALL ctl_opn( numicb, 'icebergs.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', & -1, numout, lwp, narea ) icb_alloc = 0 ALLOCATE(berg_grid, STAT=ill) icb_alloc = icb_alloc + ill ! ALLOCATE( berg_grid%calving (jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( berg_grid%calving_hflx (jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( berg_grid%stored_heat (jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( berg_grid%floating_melt(jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( berg_grid%maxclass(jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ! ALLOCATE( berg_grid%stored_ice (jpi,jpj,nclasses) , STAT=ill) icb_alloc = icb_alloc + ill ! ALLOCATE( berg_grid%tmp (jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ! ! expanded arrays for bilinear interpolation ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( vo_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( ff_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( ua_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( va_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill #if defined key_lim2 || defined key_lim3 ALLOCATE( ui_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( vi_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill #endif ALLOCATE( ssh_e(0:jpi+1,0:jpj+1) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( first_width (nclasses) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( first_length (nclasses) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( src_calving(jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( src_calving_hflx(jpi,jpj) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( nicbfldpts(jpi) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( nicbflddest(jpi) , STAT=ill) icb_alloc = icb_alloc + ill ALLOCATE( nicbfldproc(jpni) , STAT=ill) icb_alloc = icb_alloc + ill IF( lk_mpp ) CALL mpp_sum ( icb_alloc ) IF( icb_alloc > 0 ) CALL ctl_warn('icb_alloc: allocation of arrays failed') END FUNCTION icb_alloc SUBROUTINE icb_init( pdt, kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_init *** !! !! ** Purpose : iceberg initialization. !! !! ** Method : - blah blah !!---------------------------------------------------------------------- REAL(wp) , INTENT(in) :: pdt ! iceberg time-step (rdt*nn_fsbc) INTEGER , INTENT(in) :: kt ! time step number ! INTEGER :: ji, jj, jn ! loop counters INTEGER :: i1, i2, i3 ! dummy integers INTEGER :: ii, inum, ivar INTEGER :: istat1, istat2, istat3 CHARACTER(len=300) :: cl_sdist !!---------------------------------------------------------------------- CALL icb_nam ! Read and print namelist parameters IF( .NOT. ln_icebergs ) RETURN ! ! allocate gridded fields IF( icb_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' ) ! set parameters (mostly from namelist) ! berg_dt = pdt first_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) ) ) first_length(:) = rn_LoW_ratio * first_width(:) berg_grid%calving (:,:) = 0._wp berg_grid%calving_hflx (:,:) = 0._wp berg_grid%stored_heat (:,:) = 0._wp berg_grid%floating_melt(:,:) = 0._wp berg_grid%maxclass (:,:) = nclasses berg_grid%stored_ice (:,:,:) = 0._wp berg_grid%tmp (:,:) = 0._wp src_calving (:,:) = 0._wp src_calving_hflx (:,:) = 0._wp ! domain for icebergs IF( lk_mpp .AND. jpni == 1 ) & CALL ctl_stop('icbinit: having ONE processor in x currently does not work') ! for the north fold we work out which points communicate by asking ! lbc_lnk to pass processor number (valid even in single processor case) ! borrow src_calving arrays for this ! ! pack i and j together using a scaling of a power of 10 nicbpack = 10000 IF( jpiglo >= nicbpack ) CALL ctl_stop('icbini: processor index packing failure') nicbfldproc(:) = -1 DO jj = 1, jpj DO ji = 1, jpi src_calving_hflx(ji,jj) = narea src_calving(ji,jj) = nicbpack*(njmpp+jj-1) + nimpp+ji-1 ENDDO ENDDO CALL lbc_lnk( src_calving_hflx, 'T', 1._wp ) CALL lbc_lnk( src_calving , 'T', 1._wp ) ! work out interior of processor from exchange array ! first entry with narea for this processor is left hand interior index ! last entry is right hand interior index jj = jpj/2 nicbdi = -1 nicbei = -1 DO ji = 1,jpi i3 = INT( src_calving(ji,jj) ) i2 = INT( i3/nicbpack ) i1 = i3 - i2*nicbpack i3 = INT( src_calving_hflx(ji,jj) ) IF( i1 == nimpp+ji-1 .AND. i3 == narea ) THEN IF( nicbdi < 0 ) THEN nicbdi = ji ELSE nicbei = ji ENDIF ENDIF END DO ! repeat for j direction ji = jpi/2 nicbdj = -1 nicbej = -1 DO jj = 1,jpj i3 = INT( src_calving(ji,jj) ) i2 = INT( i3/nicbpack ) i1 = i3 - i2*nicbpack i3 = INT( src_calving_hflx(ji,jj) ) IF( i2 == njmpp+jj-1 .AND. i3 == narea ) THEN IF( nicbdj < 0 ) THEN nicbdj = jj ELSE nicbej = jj ENDIF ENDIF END DO ! special for east-west boundary exchange we save the destination index i1 = MAX( nicbdi-1, 1) i3 = INT( src_calving(i1,jpj/2) ) jj = INT( i3/nicbpack ) ricb_left = REAL( i3 - nicbpack*jj, wp ) i1 = MIN( nicbei+1, jpi ) i3 = INT( src_calving(i1,jpj/2) ) jj = INT( i3/nicbpack ) ricb_right = REAL( i3 - nicbpack*jj, wp ) ! north fold IF( npolj > 0 ) THEN ! ! icebergs in row nicbej+1 get passed across fold nicbfldpts(:) = INT( src_calving(:,nicbej+1) ) nicbflddest(:) = INT( src_calving_hflx(:,nicbej+1) ) ! ! work out list of unique processors to talk to DO ji = nicbdi, nicbei ii = nicbflddest(ji) DO jn = 1, jpni IF( nicbfldproc(jn) == -1 ) THEN nicbfldproc(jn) = ii EXIT ENDIF IF( nicbfldproc(jn) == ii ) EXIT ENDDO ENDDO ENDIF ! IF( nn_verbose_level > 0) THEN WRITE(numicb,*) 'processor ', narea WRITE(numicb,*) 'jpi, jpj ', jpi, jpj WRITE(numicb,*) 'nldi, nlei ', nldi, nlei WRITE(numicb,*) 'nldj, nlej ', nldj, nlej WRITE(numicb,*) 'berg i interior ', nicbdi, nicbei WRITE(numicb,*) 'berg j interior ', nicbdj, nicbej WRITE(numicb,*) 'berg left ', ricb_left WRITE(numicb,*) 'berg right ', ricb_right jj = jpj/2 WRITE(numicb,*) "central j line:" WRITE(numicb,*) "i processor" WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), ji=1,jpi) WRITE(numicb,*) "i point" WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi) ji = jpi/2 WRITE(numicb,*) "central i line:" WRITE(numicb,*) "j processor" WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), jj=1,jpj) WRITE(numicb,*) "j point" WRITE(numicb,*) (INT(src_calving(ji,jj)), jj=1,jpj) IF( npolj > 0 ) THEN WRITE(numicb,*) 'north fold destination points ' WRITE(numicb,*) nicbfldpts WRITE(numicb,*) 'north fold destination procs ' WRITE(numicb,*) nicbflddest ENDIF CALL flush(numicb) ENDIF src_calving(:,:) = 0._wp src_calving_hflx(:,:) = 0._wp ! assign each new iceberg with a unique number constructed from the processor number ! and incremented by the total number of processors num_bergs(:) = 0 num_bergs(1) = narea - jpnij ! when not generating test icebergs we need to setup calving file IF( nn_test_icebergs < 0 ) THEN ! ! maximum distribution class array does not change in time so read it once cl_sdist = TRIM( cn_dir )//TRIM( sn_icb%clname ) CALL iom_open ( cl_sdist, inum ) ! open file ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. ) IF( ivar > 0 ) THEN CALL iom_get ( inum, jpdom_data, 'maxclass', src_calving ) ! read the max distribution array berg_grid%maxclass(:,:) = INT( src_calving ) src_calving(:,:) = 0._wp ENDIF CALL iom_close( inum ) ! close file ! WRITE(numicb,*) WRITE(numicb,*) ' calving read in a file' ALLOCATE( sf_icb(1), STAT=istat1 ) ! Create sf_icb structure (calving) ALLOCATE( sf_icb(1)%fnow(jpi,jpj,1), STAT=istat2 ) ALLOCATE( sf_icb(1)%fdta(jpi,jpj,1,2), STAT=istat3 ) IF( istat1+istat2+istat3 > 0 ) THEN CALL ctl_stop( 'sbc_icb: unable to allocate sf_icb structure' ) ; RETURN ENDIF ! ! fill sf_icb with the namelist (sn_icb) and control print CALL fld_fill( sf_icb, (/ sn_icb /), cn_dir, 'icb_init', 'read calving data', 'namicb' ) ! ENDIF IF( .NOT.ln_rstart ) THEN IF( nn_test_icebergs > 0 ) CALL icb_gen() ELSE IF( nn_test_icebergs > 0 ) THEN CALL icb_gen() ELSE ! CALL icebergs_read_restart() l_restarted_bergs = .TRUE. ENDIF ENDIF ! IF( nn_sample_rate .GT. 0 ) CALL traj_init( nitend ) ! CALL icb_budget_init() ! IF( nn_verbose_level >= 2 ) CALL print_bergs('icb_init, initial status', nit000-1) ! END SUBROUTINE icb_init SUBROUTINE icb_gen() ! Local variables INTEGER :: ji, jj, ibergs TYPE(iceberg) :: localberg ! NOT a pointer but an actual local variable TYPE(point) :: localpt INTEGER :: iyr, imon, iday, ihr, imin, isec INTEGER :: iberg ! For convenience iberg = nn_test_icebergs ! call get_date(Time, iyr, imon, iday, ihr, imin, isec) ! Convert nemo time variables from dom_oce into local versions iyr = nyear imon = nmonth iday = nday ihr = INT(nsec_day/3600) imin = INT((nsec_day-ihr*3600)/60) isec = nsec_day - ihr*3600 - imin*60 ! no overlap for icebergs since we want only one instance of each across the whole domain ! so restrict area of interest ! use tmask here because tmask_i has been doctored on one side of the north fold line DO jj = nicbdj,nicbej DO ji = nicbdi,nicbei IF (tmask(ji,jj,1) .GT. 0._wp .AND. & gphit(ji,jj) .GT. rn_test_box(3) .AND. gphit(ji,jj) .LT. rn_test_box(4) .AND. & glamt(ji,jj) .GT. rn_test_box(1) .AND. glamt(ji,jj) .LT. rn_test_box(2)) THEN localberg%mass_scaling = rn_mass_scaling(iberg) localpt%xi = REAL( nimpp+ji-1, wp ) localpt%yj = REAL( njmpp+jj-1, wp ) localpt%lon = bilin(glamt, localpt%xi, localpt%yj, 'T', 0, 0 ) localpt%lat = bilin(gphit, localpt%xi, localpt%yj, 'T', 0, 0 ) localpt%mass = rn_initial_mass(iberg) localpt%thickness = rn_initial_thickness(iberg) localpt%width = first_width(iberg) localpt%length = first_length(iberg) localpt%year = iyr localpt%day = FLOAT(iday)+(FLOAT(ihr)+FLOAT(imin)/60._wp)/24._wp localpt%mass_of_bits = 0._wp localpt%heat_density = 0._wp localpt%uvel = 0._wp localpt%vvel = 0._wp CALL increment_kounter() localberg%number(:) = num_bergs(:) call add_new_berg_to_list(localberg, localpt) ENDIF ENDDO ENDDO ibergs = count_bergs() IF( lk_mpp ) CALL mpp_sum(ibergs) WRITE(numicb,'(a,i6,a)') 'diamonds, icb_gen: ',ibergs,' were generated' END SUBROUTINE icb_gen SUBROUTINE icb_nam !!---------------------------------------------------------------------- !! *** ROUTINE icb_nam *** !! !! ** Purpose : read domaine namelists and print the variables. !! !! ** input : - namberg namelist !!---------------------------------------------------------------------- NAMELIST/namberg/ ln_icebergs, ln_bergdia, nn_sample_rate, rn_initial_mass, & & rn_distribution, rn_mass_scaling, rn_initial_thickness, nn_verbose_write, & & rn_rho_bergs, rn_LoW_ratio, nn_verbose_level, ln_operator_splitting, & & rn_bits_erosion_fraction, rn_sicn_shift, ln_passive_mode, & & ln_time_average_weight, nn_test_icebergs, rn_test_box, rn_speed_limit, & & cn_dir, sn_icb INTEGER :: jn ! loop counter REAL(wp) :: zfact !!---------------------------------------------------------------------- ! (NB: frequency positive => hours, negative => months) ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! sn_icb = FLD_N( 'calving' , -1 , 'calving' , .TRUE. , .TRUE. , 'yearly' , '' , '' ) REWIND( numnam ) ! Namelist namrun : iceberg parameters READ ( numnam, namberg ) IF( .NOT. ln_icebergs .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read' WRITE(numout,*) '~~~~~~~~~~~~~ ' WRITE(numout,*) 'NO icebergs used' RETURN ENDIF IF( nn_test_icebergs > nclasses ) THEN IF( lwp ) WRITE(numout,*) 'Resetting nn_test_icebergs to ',nclasses nn_test_icebergs = nclasses ENDIF zfact = SUM( rn_distribution ) IF( zfact < 1._wp ) THEN IF( zfact <= 0._wp ) THEN CALL ctl_stop('icb_init: sum of berg distribution equal to zero') ELSE rn_distribution(:) = rn_distribution(:) / zfact CALL ctl_warn('icb_init: sum of berg input distribution not equal to one and so RESCALED') ENDIF ENDIF IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read' WRITE(numout,*) '~~~~~~~~~~~~~ ' WRITE(numout,*) ' Calculate budgets ln_bergdia = ', ln_bergdia WRITE(numout,*) ' Period between sampling of position for trajectory storage nn_sample_rate = ', nn_sample_rate WRITE(numout,*) ' Mass thresholds between iceberg classes (kg) rn_initial_mass =' DO jn=1,nclasses WRITE(numout,'(a,f15.2)') ' ',rn_initial_mass(jn) ENDDO WRITE(numout,*) ' Fraction of calving to apply to this class (non-dim) rn_distribution =' DO jn=1,nclasses WRITE(numout,'(a,f10.2)') ' ',rn_distribution(jn) ENDDO WRITE(numout,*) ' Ratio between effective and real iceberg mass (non-dim) rn_mass_scaling = ' DO jn=1,nclasses WRITE(numout,'(a,f10.2)') ' ',rn_mass_scaling(jn) ENDDO WRITE(numout,*) ' Total thickness of newly calved bergs (m) rn_initial_thickness = ' DO jn=1,nclasses WRITE(numout,'(a,f10.2)') ' ',rn_initial_thickness(jn) ENDDO WRITE(numout,*) ' Timesteps between verbose messages nn_verbose_write = ', nn_verbose_write WRITE(numout,*) ' Density of icebergs rn_rho_bergs = ', rn_rho_bergs WRITE(numout,*) ' Initial ratio L/W for newly calved icebergs rn_LoW_ratio = ', rn_LoW_ratio WRITE(numout,*) ' Turn on more verbose output level = ', nn_verbose_level WRITE(numout,*) ' Use first order operator splitting for thermodynamics ', & & 'use_operator_splitting = ', ln_operator_splitting WRITE(numout,*) ' Fraction of erosion melt flux to divert to bergy bits ', & & 'bits_erosion_fraction = ', rn_bits_erosion_fraction WRITE(numout,*) ' Shift of sea-ice concentration in erosion flux modulation ', & & '(0