LCOV - code coverage report
Current view: top level - configuration/driver - icedrv_init.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 553 639 86.54 %
Date: 2020-07-03 20:11:40 Functions: 5 5 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! parameter and variable initializations
       4             : !
       5             : ! authors Elizabeth C. Hunke, LANL
       6             : 
       7             :       module icedrv_init
       8             : 
       9             :       use icedrv_kinds
      10             :       use icedrv_constants, only: nu_diag, ice_stdout, nu_diag_out, nu_nml
      11             :       use icedrv_constants, only: c0, c1, c2, c3, p2, p5
      12             :       use icedrv_domain_size, only: nx
      13             :       use icepack_intfc, only: icepack_init_parameters
      14             :       use icepack_intfc, only: icepack_init_fsd
      15             :       use icepack_intfc, only: icepack_init_tracer_flags
      16             :       use icepack_intfc, only: icepack_init_tracer_sizes
      17             :       use icepack_intfc, only: icepack_init_tracer_indices
      18             :       use icepack_intfc, only: icepack_init_trcr
      19             :       use icepack_intfc, only: icepack_query_parameters
      20             :       use icepack_intfc, only: icepack_query_tracer_flags
      21             :       use icepack_intfc, only: icepack_query_tracer_sizes
      22             :       use icepack_intfc, only: icepack_query_tracer_indices
      23             :       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
      24             :       use icedrv_system, only: icedrv_system_abort
      25             : 
      26             :       implicit none
      27             :       private
      28             :       public :: input_data, init_grid2, init_state, init_fsd
      29             : 
      30             :       character(len=char_len_long), public :: &
      31             :          ice_ic      ! method of ice cover initialization
      32             :                      ! 'default' or 'none' => conditions specified in code
      33             :                      ! restart = .true. overwrites default initial
      34             :                      !    condition using filename given by ice_ic
      35             : 
      36             :       real (kind=dbl_kind), dimension (nx), public :: &
      37             :          TLON   , & ! longitude of temp pts (radians)
      38             :          TLAT       ! latitude of temp pts (radians)
      39             : 
      40             :       logical (kind=log_kind), &
      41             :          dimension (nx), public :: &
      42             :          tmask  , & ! land/boundary mask, thickness (T-cell)
      43             :          lmask_n, & ! northern hemisphere mask
      44             :          lmask_s    ! southern hemisphere mask
      45             : 
      46             : !=======================================================================
      47             : 
      48             :       contains
      49             : 
      50             : !=======================================================================
      51             : 
      52             : ! Namelist variables, set to default values; may be altered
      53             : ! at run time
      54             : !
      55             : ! author Elizabeth C. Hunke, LANL
      56             : 
      57          45 :       subroutine input_data
      58             : 
      59             :       use icedrv_diagnostics, only: diag_file, nx_names
      60             :       use icedrv_domain_size, only: nilyr, nslyr, nblyr, max_ntrcr, ncat
      61             :       use icedrv_domain_size, only: n_iso, n_aero, nfsd
      62             :       use icedrv_calendar, only: year_init, istep0
      63             :       use icedrv_calendar, only: dumpfreq, diagfreq, dump_last
      64             :       use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years
      65             :       use icedrv_restart_shared, only: restart, restart_dir, restart_file
      66             :       use icedrv_flux, only: update_ocn_f, l_mpond_fresh, cpl_bgc
      67             :       use icedrv_flux, only: default_season
      68             :       use icedrv_forcing, only: precip_units,    fyear_init,      ycycle
      69             :       use icedrv_forcing, only: atm_data_type,   ocn_data_type,   bgc_data_type
      70             :       use icedrv_forcing, only: atm_data_file,   ocn_data_file,   bgc_data_file
      71             :       use icedrv_forcing, only: ice_data_file
      72             :       use icedrv_forcing, only: atm_data_format, ocn_data_format, bgc_data_format
      73             :       use icedrv_forcing, only: data_dir
      74             :       use icedrv_forcing, only: oceanmixed_ice, restore_ocn, trestore
      75             : 
      76             :       ! local variables
      77             : 
      78             :       character (32) :: &
      79             :          nml_filename = 'icepack_in' ! namelist input file name
      80             : 
      81             :       integer (kind=int_kind) :: &
      82             :          nml_error, & ! namelist i/o error flag
      83             :          n            ! loop index
      84             : 
      85             :       character (len=char_len) :: diag_file_names
      86             : 
      87          32 :       real (kind=dbl_kind) :: ustar_min, albicev, albicei, albsnowv, albsnowi, &
      88          80 :          ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, &
      89          80 :          mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
      90          64 :          a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
      91          48 :          phi_c_slow_mode, phi_i_mushy, kalg, emissivity
      92             : 
      93             :       integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, &
      94             :          natmiter, kitd, kcatbound
      95             : 
      96             :       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
      97             :          tfrz_option, frzpnd, atmbndy, wave_spec_type
      98             : 
      99             :       ! Flux convergence tolerance
     100          16 :       real (kind=dbl_kind) :: atmiter_conv
     101             : 
     102             :       logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair
     103             :       logical (kind=log_kind) :: conserv_check
     104             : 
     105             :       integer (kind=int_kind) :: ntrcr
     106             :       logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond
     107             :       logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd
     108             :       logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo, wave_spec
     109             :       integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY
     110             :       integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, &
     111             :                                  nt_aero, nt_fsd, nt_isosno, nt_isoice
     112             : 
     113          16 :       real (kind=real_kind) :: rpcesm, rplvl, rptopo 
     114          32 :       real (kind=dbl_kind) :: Cf, puny
     115             :       character(len=*), parameter :: subname='(input_data)'
     116             : 
     117             :       !-----------------------------------------------------------------
     118             :       ! Namelist variables
     119             :       !-----------------------------------------------------------------
     120             : 
     121             :       namelist /setup_nml/ &
     122             :         days_per_year,  use_leap_years, year_init,       istep0,        &
     123             :         dt,             npt,            ndtd,            dump_last,     &
     124             :         ice_ic,         restart,        restart_dir,     restart_file,  &
     125             :         dumpfreq,       diagfreq,       diag_file,       cpl_bgc,       &
     126             :         conserv_check
     127             : 
     128             :       namelist /grid_nml/ &
     129             :         kcatbound
     130             : 
     131             :       namelist /thermo_nml/ &
     132             :         kitd,           ktherm,          ksno,     conduct,             &
     133             :         a_rapid_mode,   Rac_rapid_mode,  aspect_rapid_mode,             &
     134             :         dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy
     135             : 
     136             :       namelist /dynamics_nml/ &
     137             :         kstrength,      krdg_partic,    krdg_redist,    mu_rdg,         &
     138             :         Cf
     139             : 
     140             :       namelist /shortwave_nml/ &
     141             :         shortwave,      albedo_type,                                    &
     142             :         albicev,        albicei,         albsnowv,      albsnowi,       &
     143             :         ahmax,          R_ice,           R_pnd,         R_snw,          &
     144             :         dT_mlt,         rsnw_mlt,        kalg
     145             : 
     146             :       namelist /ponds_nml/ &
     147             :         hs0,            dpscale,         frzpnd,                        &
     148             :         rfracmin,       rfracmax,        pndaspect,     hs1,            &
     149             :         hp1
     150             : 
     151             :       namelist /forcing_nml/ &
     152             :         atmbndy,         calc_strair,     calc_Tsfc,       &
     153             :         update_ocn_f,    l_mpond_fresh,   ustar_min,       &
     154             :         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
     155             :         formdrag,        highfreq,        natmiter,        &
     156             :         atmiter_conv,                        &
     157             :         tfrz_option,     default_season,  wave_spec_type,  &
     158             :         precip_units,    fyear_init,      ycycle,          &
     159             :         atm_data_type,   ocn_data_type,   bgc_data_type,   &
     160             :         atm_data_file,   ocn_data_file,   bgc_data_file,   &
     161             :         ice_data_file,                                     &
     162             :         atm_data_format, ocn_data_format, bgc_data_format, &
     163             :         data_dir,        trestore,        restore_ocn
     164             : 
     165             :       namelist /tracer_nml/   &
     166             :         tr_iage,      &
     167             :         tr_FY,        &
     168             :         tr_lvl,       &
     169             :         tr_pond_cesm, &
     170             :         tr_pond_lvl,  &
     171             :         tr_pond_topo, &
     172             :         tr_aero,      &
     173             :         tr_fsd,       &
     174             :         tr_iso
     175             : 
     176             :       !-----------------------------------------------------------------
     177             :       ! query Icepack values
     178             :       !-----------------------------------------------------------------
     179             : 
     180             :       call icepack_query_parameters(ustar_min_out=ustar_min, Cf_out=Cf, &
     181             :            albicev_out=albicev, albicei_out=albicei, ksno_out = ksno,   &
     182             :            albsnowv_out=albsnowv, albsnowi_out=albsnowi, &
     183             :            natmiter_out=natmiter, ahmax_out=ahmax, shortwave_out=shortwave, &
     184             :            atmiter_conv_out = atmiter_conv, &
     185             :            albedo_type_out=albedo_type, R_ice_out=R_ice, R_pnd_out=R_pnd, &
     186             :            R_snw_out=R_snw, dT_mlt_out=dT_mlt, rsnw_mlt_out=rsnw_mlt, &
     187             :            kstrength_out=kstrength, krdg_partic_out=krdg_partic, &
     188             :            krdg_redist_out=krdg_redist, mu_rdg_out=mu_rdg, &
     189             :            atmbndy_out=atmbndy, calc_strair_out=calc_strair, &
     190             :            formdrag_out=formdrag, highfreq_out=highfreq, &
     191             :            emissivity_out=emissivity, &
     192             :            kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, & 
     193             :            dpscale_out=dpscale, frzpnd_out=frzpnd, &
     194             :            rfracmin_out=rfracmin, rfracmax_out=rfracmax, &
     195             :            pndaspect_out=pndaspect, hs1_out=hs1, hp1_out=hp1, &
     196             :            ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc, &
     197             :            update_ocn_f_out = update_ocn_f, &
     198             :            conduct_out=conduct, a_rapid_mode_out=a_rapid_mode, &
     199             :            Rac_rapid_mode_out=Rac_rapid_mode, &
     200             :            aspect_rapid_mode_out=aspect_rapid_mode, &
     201             :            dSdt_slow_mode_out=dSdt_slow_mode, &
     202             :            phi_c_slow_mode_out=phi_c_slow_mode, &
     203             :            phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
     204             :            tfrz_option_out=tfrz_option, kalg_out=kalg, &
     205             :            fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
     206          45 :            wave_spec_type_out=wave_spec_type)
     207          45 :       call icepack_warnings_flush(nu_diag)
     208          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     209           0 :           file=__FILE__, line=__LINE__)
     210             : 
     211             :       !-----------------------------------------------------------------
     212             :       ! default values
     213             :       !-----------------------------------------------------------------
     214             : 
     215          45 :       days_per_year = 365    ! number of days in a year
     216          45 :       use_leap_years= .false.! if true, use leap years (Feb 29)
     217          45 :       year_init = 0          ! initial year
     218          45 :       istep0 = 0             ! no. of steps taken in previous integrations,
     219             :                              ! real (dumped) or imagined (to set calendar)
     220          45 :       dt = 3600.0_dbl_kind   ! time step, s      
     221          45 :       npt = 99999            ! total number of time steps (dt) 
     222          45 :       diagfreq = 24          ! how often diag output is written
     223          45 :       diag_file = 'ice_diag' ! history file name prefix
     224          45 :       cpl_bgc = .false.      ! 
     225          45 :       dumpfreq='y'           ! restart frequency option
     226          45 :       dump_last=.false.      ! restart at end of run
     227          45 :       restart = .false.      ! if true, read restart files for initialization
     228          45 :       restart_dir  = './'    ! write to executable dir for default
     229          45 :       restart_file = 'iced'  ! restart file name prefix
     230          45 :       ice_ic       = 'default'      ! initial conditions are specified in the code
     231             :                                     ! otherwise, the filename for reading restarts
     232          45 :       ndtd = 1               ! dynamic time steps per thermodynamic time step
     233          45 :       l_mpond_fresh = .false.     ! logical switch for including meltpond freshwater
     234             :                                   ! flux feedback to ocean model
     235          45 :       default_season  = 'winter'  ! default forcing data, if data is not read in
     236          45 :       fyear_init      = 1998      ! initial forcing year
     237          45 :       ycycle          = 1         ! number of years in forcing cycle
     238          45 :       atm_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
     239          45 :       atm_data_type   = 'default' ! source of atmospheric forcing data
     240          45 :       atm_data_file   = ' '       ! atmospheric forcing data file
     241          45 :       precip_units    = 'mks'     ! 'mm_per_month' or
     242             :                                   ! 'mm_per_sec' = 'mks' = kg/m^2 s
     243          45 :       oceanmixed_ice  = .false.   ! if true, use internal ocean mixed layer
     244          45 :       wave_spec_type  = 'none'    ! type of wave spectrum forcing
     245          45 :       ocn_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
     246          45 :       ocn_data_type   = 'default' ! source of ocean forcing data
     247          45 :       ocn_data_file   = ' '       ! ocean forcing data file
     248          45 :       ice_data_file   = ' '       ! ice forcing data file (opening, closing)
     249          45 :       bgc_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
     250          45 :       bgc_data_type   = 'default' ! source of BGC forcing data
     251          45 :       bgc_data_file   = ' '       ! biogeochemistry forcing data file
     252          45 :       data_dir    = ' '           ! root location of data files
     253          45 :       restore_ocn     = .false.   ! restore sst if true
     254          45 :       trestore        = 90        ! restoring timescale, days (0 instantaneous)
     255             : 
     256             :       ! extra tracers
     257          45 :       tr_iage      = .false. ! ice age
     258          45 :       tr_FY        = .false. ! ice age
     259          45 :       tr_lvl       = .false. ! level ice 
     260          45 :       tr_pond_cesm = .false. ! CESM melt ponds
     261          45 :       tr_pond_lvl  = .false. ! level-ice melt ponds
     262          45 :       tr_pond_topo = .false. ! explicit melt ponds (topographic)
     263          45 :       tr_aero      = .false. ! aerosols
     264          45 :       tr_fsd       = .false. ! floe size distribution
     265          45 :       tr_iso       = .false. ! isotopes
     266             : 
     267             :       !-----------------------------------------------------------------
     268             :       ! read from input file
     269             :       !-----------------------------------------------------------------
     270             : 
     271          45 :       open (nu_nml, file=nml_filename, status='old',iostat=nml_error)
     272          45 :       if (nml_error /= 0) then
     273           0 :          nml_error = -1
     274             :       else
     275          45 :          nml_error =  1
     276             :       endif
     277             :       
     278          90 :       do while (nml_error > 0)
     279          45 :          print*,'Reading namelist file   ',nml_filename
     280             : 
     281          45 :          print*,'Reading setup_nml'
     282          45 :          read(nu_nml, nml=setup_nml,iostat=nml_error)
     283          45 :          if (nml_error /= 0) exit
     284             : 
     285          45 :          print*,'Reading grid_nml'
     286          45 :          read(nu_nml, nml=grid_nml,iostat=nml_error)
     287          45 :          if (nml_error /= 0) exit
     288             : 
     289          45 :          print*,'Reading tracer_nml'
     290          45 :          read(nu_nml, nml=tracer_nml,iostat=nml_error)
     291          45 :          if (nml_error /= 0) exit
     292             : 
     293          45 :          print*,'Reading thermo_nml'
     294          45 :          read(nu_nml, nml=thermo_nml,iostat=nml_error)
     295          45 :          if (nml_error /= 0) exit
     296             : 
     297          45 :          print*,'Reading shortwave_nml'
     298          45 :          read(nu_nml, nml=shortwave_nml,iostat=nml_error)
     299          45 :          if (nml_error /= 0) exit
     300             : 
     301          45 :          print*,'Reading ponds_nml'
     302          45 :          read(nu_nml, nml=ponds_nml,iostat=nml_error)
     303          45 :          if (nml_error /= 0) exit
     304             : 
     305          45 :          print*,'Reading forcing_nml'
     306          45 :          read(nu_nml, nml=forcing_nml,iostat=nml_error)
     307          90 :          if (nml_error /= 0) exit
     308             :       end do
     309          45 :       if (nml_error == 0) close(nu_nml)
     310          45 :       if (nml_error /= 0) then
     311           0 :          write(ice_stdout,*) 'error reading namelist'
     312           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     313             :       endif
     314          45 :       close(nu_nml)
     315             :       
     316             :       !-----------------------------------------------------------------
     317             :       ! set up diagnostics output and resolve conflicts
     318             :       !-----------------------------------------------------------------
     319             :       
     320          45 :       write(ice_stdout,*) 'Diagnostic output will be in files '
     321          45 :       write(ice_stdout,*)'    ','icepack.runlog.timestamp'
     322             : 
     323         225 :       do n = 1,nx
     324         225 :          write(nx_names(n),'(a,i2.2)') 'point_',n
     325             :       enddo      
     326          45 :       nx_names(1) = 'icefree'
     327          45 :       nx_names(2) = 'slab'
     328          45 :       nx_names(3) = 'full_ITD'
     329          45 :       nx_names(4) = 'land'
     330             : 
     331         225 :       do n = 1,nx
     332         180 :          diag_file_names=' '
     333         180 :          write(diag_file_names,'(a,a,a)') trim(diag_file),'.',trim(nx_names(n))
     334         180 :          write(ice_stdout,*)'    ',trim(diag_file_names)
     335         225 :          open(nu_diag_out+n-1, file=diag_file_names, status='unknown')
     336             :       end do
     337             :       
     338          45 :       write(nu_diag,*) '-----------------------------------'
     339          45 :       write(nu_diag,*) '  ICEPACK model diagnostic output  '
     340          45 :       write(nu_diag,*) '-----------------------------------'
     341          45 :       write(nu_diag,*) ' '
     342             : 
     343           3 :       if (ncat == 1 .and. kitd == 1) then
     344           0 :          write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1.'
     345           0 :          write (nu_diag,*) 'Use kitd = 0 (delta function ITD) with kcatbound = 0'
     346           0 :          write (nu_diag,*) 'or for column configurations use kcatbound = -1'
     347           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     348             :       endif
     349             : 
     350          42 :       if (ncat /= 1 .and. kcatbound == -1) then
     351           0 :          write (nu_diag,*) 'WARNING: ITD required for ncat > 1'
     352           0 :          write (nu_diag,*) 'WARNING: Setting kitd and kcatbound to default values'
     353           0 :          kitd = 1
     354           0 :          kcatbound = 0
     355             :       endif
     356             : 
     357          45 :       rpcesm = c0
     358          45 :       rplvl  = c0
     359          45 :       rptopo = c0
     360          45 :       if (tr_pond_cesm) rpcesm = c1
     361          45 :       if (tr_pond_lvl ) rplvl  = c1
     362          45 :       if (tr_pond_topo) rptopo = c1
     363             : 
     364          45 :       tr_pond = .false. ! explicit melt ponds
     365          45 :       if (rpcesm + rplvl + rptopo > puny) tr_pond = .true.
     366             : 
     367          45 :       if (rpcesm + rplvl + rptopo > c1 + puny) then
     368           0 :          write (nu_diag,*) 'WARNING: Must use only one melt pond scheme'
     369           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     370             :       endif
     371             : 
     372          45 :       if (tr_pond_lvl .and. .not. tr_lvl) then
     373           0 :          write (nu_diag,*) 'WARNING: tr_pond_lvl=T but tr_lvl=F'
     374           0 :          write (nu_diag,*) 'WARNING: Setting tr_lvl=T'
     375           0 :          tr_lvl = .true.
     376             :       endif
     377             : 
     378          45 :       if (tr_pond_lvl .and. abs(hs0) > puny) then
     379           0 :          write (nu_diag,*) 'WARNING: tr_pond_lvl=T and hs0/=0'
     380           0 :          write (nu_diag,*) 'WARNING: Setting hs0=0'
     381           0 :          hs0 = c0
     382             :       endif
     383             : 
     384          45 :       if (tr_pond_cesm .and. trim(frzpnd) /= 'cesm') then
     385           2 :          write (nu_diag,*) 'WARNING: tr_pond_cesm=T'
     386           2 :          write (nu_diag,*) 'WARNING: frzpnd, dpscale not used'
     387           2 :          frzpnd = 'cesm'
     388             :       endif
     389             : 
     390          45 :       if (trim(shortwave) /= 'dEdd' .and. tr_pond .and. calc_tsfc) then
     391           0 :          write (nu_diag,*) 'WARNING: Must use dEdd shortwave'
     392           0 :          write (nu_diag,*) 'WARNING: with tr_pond and calc_tsfc=T.'
     393           0 :          write (nu_diag,*) 'WARNING: Setting shortwave = dEdd'
     394           0 :          shortwave = 'dEdd'
     395             :       endif
     396             : 
     397          43 :       if (tr_iso .and. n_iso==0) then
     398           0 :          write (nu_diag,*) 'WARNING: isotopes activated but'
     399           0 :          write (nu_diag,*) 'WARNING: not allocated in tracer array.'
     400           0 :          write (nu_diag,*) 'WARNING: Activate in compilation script.'
     401           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     402             :       endif
     403             : 
     404           5 :       if (tr_aero .and. n_aero==0) then
     405           0 :          write (nu_diag,*) 'WARNING: aerosols activated but'
     406           0 :          write (nu_diag,*) 'WARNING: not allocated in tracer array.'
     407           0 :          write (nu_diag,*) 'WARNING: Activate in compilation script.'
     408           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     409             :       endif
     410             : 
     411          45 :       if (tr_aero .and. trim(shortwave) /= 'dEdd') then
     412           0 :          write (nu_diag,*) 'WARNING: aerosols activated but dEdd'
     413           0 :          write (nu_diag,*) 'WARNING: shortwave is not.'
     414           0 :          write (nu_diag,*) 'WARNING: Setting shortwave = dEdd'
     415           0 :          shortwave = 'dEdd'
     416             :       endif
     417             : 
     418          45 :       rfracmin = min(max(rfracmin,c0),c1)
     419          45 :       rfracmax = min(max(rfracmax,c0),c1)
     420             : 
     421          45 :       if (ktherm == 2 .and. .not. calc_Tsfc) then
     422           0 :          write (nu_diag,*) 'WARNING: ktherm = 2 and calc_Tsfc = F'
     423           0 :          write (nu_diag,*) 'WARNING: Setting calc_Tsfc = T'
     424           0 :          calc_Tsfc = .true.
     425             :       endif
     426             : 
     427          45 :       if (ktherm == 1 .and. trim(tfrz_option) /= 'linear_salt') then
     428             :          write (nu_diag,*) &
     429           0 :          'WARNING: ktherm = 1 and tfrz_option = ',trim(tfrz_option)
     430             :          write (nu_diag,*) &
     431           0 :          'WARNING: For consistency, set tfrz_option = linear_salt'
     432             :       endif
     433          45 :       if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then
     434             :          write (nu_diag,*) &
     435           0 :          'WARNING: ktherm = 2 and tfrz_option = ',trim(tfrz_option)
     436             :          write (nu_diag,*) &
     437           0 :          'WARNING: For consistency, set tfrz_option = mushy'
     438             :       endif
     439             : 
     440          45 :       if (formdrag) then
     441           3 :       if (trim(atmbndy) == 'constant') then
     442           0 :          write (nu_diag,*) 'WARNING: atmbndy = constant not allowed with formdrag'
     443           0 :          write (nu_diag,*) 'WARNING: Setting atmbndy = default'
     444           0 :          atmbndy = 'default'
     445             :       endif
     446             : 
     447           3 :       if (.not. calc_strair) then
     448           0 :          write (nu_diag,*) 'WARNING: formdrag=T but calc_strair=F'
     449           0 :          write (nu_diag,*) 'WARNING: Setting calc_strair=T'
     450           0 :          calc_strair = .true.
     451             :       endif
     452             : 
     453           3 :       if (tr_pond_cesm) then
     454           0 :          write (nu_diag,*) 'ERROR: formdrag=T but frzpnd=cesm' 
     455           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
     456             :       endif
     457             : 
     458           3 :       if (.not. tr_lvl) then
     459           3 :          write (nu_diag,*) 'WARNING: formdrag=T but tr_lvl=F'
     460           3 :          write (nu_diag,*) 'WARNING: Setting tr_lvl=T'
     461           3 :          tr_lvl = .true.
     462             :       endif
     463             :       endif
     464             : 
     465          45 :       if (trim(fbot_xfer_type) == 'Cdn_ocn' .and. .not. formdrag)  then
     466           0 :          write (nu_diag,*) 'WARNING: formdrag=F but fbot_xfer_type=Cdn_ocn'
     467           0 :          write (nu_diag,*) 'WARNING: Setting fbot_xfer_type = constant'
     468           0 :          fbot_xfer_type = 'constant'
     469             :       endif
     470             : 
     471          45 :       wave_spec = .false.
     472          45 :       if (tr_fsd .and. (trim(wave_spec_type) /= 'none')) wave_spec = .true.
     473             : 
     474             :       !-----------------------------------------------------------------
     475             :       ! spew
     476             :       !-----------------------------------------------------------------
     477             : 
     478          45 :          write(nu_diag,*) ' Document ice_in namelist parameters:'
     479          45 :          write(nu_diag,*) ' ==================================== '
     480          45 :          write(nu_diag,*) ' '
     481          45 :          write(nu_diag,1020) ' days_per_year             = ', days_per_year
     482          45 :          write(nu_diag,1010) ' use_leap_years            = ', use_leap_years
     483          45 :          write(nu_diag,1020) ' year_init                 = ', year_init
     484          45 :          write(nu_diag,1020) ' istep0                    = ', istep0
     485          45 :          write(nu_diag,1000) ' dt                        = ', dt
     486          45 :          write(nu_diag,1020) ' npt                       = ', npt
     487          45 :          write(nu_diag,1020) ' diagfreq                  = ', diagfreq
     488          45 :          write(nu_diag,1030) ' dumpfreq                  = ', &
     489          90 :                                trim(dumpfreq)
     490          45 :          write(nu_diag,1010) ' dump_last                 = ', dump_last
     491          45 :          write(nu_diag,1010) ' restart                   = ', restart
     492          45 :          write(nu_diag,*)    ' restart_dir               = ', &
     493          90 :                                trim(restart_dir)
     494          45 :          write(nu_diag,*)    ' restart_file              = ', &
     495          90 :                                trim(restart_file)
     496          45 :          write(nu_diag,*)    ' ice_ic                    = ', &
     497          90 :                                trim(ice_ic)
     498          45 :          write(nu_diag,1010) ' conserv_check             = ', conserv_check
     499          45 :          write(nu_diag,1020) ' kitd                      = ', kitd
     500          45 :          write(nu_diag,1020) ' kcatbound                 = ', &
     501          90 :                                kcatbound
     502          45 :          write(nu_diag,1020) ' ndtd                      = ', ndtd
     503          45 :          write(nu_diag,1020) ' kstrength                 = ', kstrength
     504          45 :          write(nu_diag,1020) ' krdg_partic               = ', &
     505          90 :                                krdg_partic
     506          45 :          write(nu_diag,1020) ' krdg_redist               = ', &
     507          90 :                                krdg_redist
     508          45 :          if (krdg_redist == 1) &
     509          45 :          write(nu_diag,1000) ' mu_rdg                    = ', mu_rdg
     510          45 :          if (kstrength == 1) &
     511          45 :          write(nu_diag,1000) ' Cf                        = ', Cf
     512          45 :          write(nu_diag,1000) ' ksno                      = ', ksno
     513          45 :          write(nu_diag,1030) ' shortwave                 = ', &
     514          90 :                                trim(shortwave)
     515          45 :          if (cpl_bgc) then
     516           0 :              write(nu_diag,1000) ' BGC coupling is switched ON'
     517             :          else
     518          45 :              write(nu_diag,1000) ' BGC coupling is switched OFF'
     519             :          endif
     520             : 
     521          45 :          if (trim(shortwave) == 'dEdd') then
     522          39 :          write(nu_diag,1000) ' R_ice                     = ', R_ice
     523          39 :          write(nu_diag,1000) ' R_pnd                     = ', R_pnd
     524          39 :          write(nu_diag,1000) ' R_snw                     = ', R_snw
     525          39 :          write(nu_diag,1000) ' dT_mlt                    = ', dT_mlt
     526          39 :          write(nu_diag,1000) ' rsnw_mlt                  = ', rsnw_mlt
     527          39 :          write(nu_diag,1000) ' kalg                      = ', kalg
     528          39 :          write(nu_diag,1000) ' hp1                       = ', hp1
     529          39 :          write(nu_diag,1000) ' hs0                       = ', hs0
     530             :          else
     531           6 :          write(nu_diag,1030) ' albedo_type               = ', &
     532          12 :                                trim(albedo_type)
     533           6 :          write(nu_diag,1000) ' albicev                   = ', albicev
     534           6 :          write(nu_diag,1000) ' albicei                   = ', albicei
     535           6 :          write(nu_diag,1000) ' albsnowv                  = ', albsnowv
     536           6 :          write(nu_diag,1000) ' albsnowi                  = ', albsnowi
     537           6 :          write(nu_diag,1000) ' ahmax                     = ', ahmax
     538             :          endif
     539             : 
     540          45 :          write(nu_diag,1000) ' rfracmin                  = ', rfracmin
     541          45 :          write(nu_diag,1000) ' rfracmax                  = ', rfracmax
     542          45 :          if (tr_pond_lvl) then
     543          29 :          write(nu_diag,1000) ' hs1                       = ', hs1
     544          29 :          write(nu_diag,1000) ' dpscale                   = ', dpscale
     545          29 :          write(nu_diag,1030) ' frzpnd                    = ', trim(frzpnd)
     546             :          endif
     547          45 :          if (tr_pond .and. .not. tr_pond_lvl) &
     548           7 :          write(nu_diag,1000) ' pndaspect                 = ', pndaspect
     549             : 
     550          45 :          write(nu_diag,1020) ' ktherm                    = ', ktherm
     551          45 :          if (ktherm == 1) &
     552           9 :          write(nu_diag,1030) ' conduct                   = ', conduct
     553          45 :          write(nu_diag,1005) ' emissivity                = ', emissivity
     554          45 :          if (ktherm == 2) then
     555          30 :          write(nu_diag,1005) ' a_rapid_mode              = ', a_rapid_mode
     556          30 :          write(nu_diag,1005) ' Rac_rapid_mode            = ', Rac_rapid_mode
     557          30 :          write(nu_diag,1005) ' aspect_rapid_mode         = ', aspect_rapid_mode
     558          30 :          write(nu_diag,1005) ' dSdt_slow_mode            = ', dSdt_slow_mode
     559          30 :          write(nu_diag,1005) ' phi_c_slow_mode           = ', phi_c_slow_mode
     560          30 :          write(nu_diag,1005) ' phi_i_mushy               = ', phi_i_mushy
     561             :          endif
     562             : 
     563          45 :          write(nu_diag,1030) ' atmbndy                   = ', &
     564          90 :                                trim(atmbndy)
     565          45 :          write(nu_diag,1010) ' formdrag                  = ', formdrag
     566          45 :          write(nu_diag,1010) ' highfreq                  = ', highfreq
     567          45 :          write(nu_diag,1020) ' natmiter                  = ', natmiter
     568          45 :          write(nu_diag,1005) ' atmiter_conv              = ', atmiter_conv
     569          45 :          write(nu_diag,1010) ' calc_strair               = ', calc_strair
     570          45 :          write(nu_diag,1010) ' calc_Tsfc                 = ', calc_Tsfc
     571             : 
     572          45 :          write(nu_diag,*)    ' atm_data_type             = ', &
     573          90 :                                trim(atm_data_type)
     574          45 :          write(nu_diag,*)    ' ocn_data_type             = ', &
     575          90 :                                trim(ocn_data_type)
     576          45 :          write(nu_diag,*)    ' bgc_data_type             = ', &
     577          90 :                                trim(bgc_data_type)
     578             : 
     579          45 :          write(nu_diag,*)    ' atm_data_file             = ', &
     580          90 :                                trim(atm_data_file)
     581          45 :          write(nu_diag,*)    ' ocn_data_file             = ', &
     582          90 :                                trim(ocn_data_file)
     583          45 :          write(nu_diag,*)    ' bgc_data_file             = ', &
     584          90 :                                trim(bgc_data_file)
     585          45 :          write(nu_diag,*)    ' ice_data_file             = ', &
     586          90 :                                trim(ice_data_file)
     587             : 
     588          45 :          if (trim(atm_data_type)=='default') &
     589           9 :          write(nu_diag,*)    ' default_season            = ', trim(default_season)
     590             : 
     591          45 :          write(nu_diag,1010) ' update_ocn_f              = ', update_ocn_f
     592          45 :          write(nu_diag,1010) ' wave_spec                 = ', wave_spec
     593          45 :          if (wave_spec) &
     594           3 :          write(nu_diag,*)    ' wave_spec_type            = ', wave_spec_type
     595          45 :          write(nu_diag,1010) ' l_mpond_fresh             = ', l_mpond_fresh
     596          45 :          write(nu_diag,1005) ' ustar_min                 = ', ustar_min
     597          45 :          write(nu_diag,*)    ' fbot_xfer_type            = ', &
     598          90 :                                trim(fbot_xfer_type)
     599          45 :          write(nu_diag,1010) ' oceanmixed_ice            = ', oceanmixed_ice
     600          45 :          write(nu_diag,*)    ' tfrz_option               = ', &
     601          90 :                                trim(tfrz_option)
     602          45 :          write(nu_diag,1010) ' restore_ocn               = ', restore_ocn
     603          45 :          if (restore_ocn) &
     604          11 :          write(nu_diag,1005) ' trestore                  = ', trestore
     605             : 
     606             :          ! tracers
     607          45 :          write(nu_diag,1010) ' tr_iage                   = ', tr_iage
     608          45 :          write(nu_diag,1010) ' tr_FY                     = ', tr_FY
     609          45 :          write(nu_diag,1010) ' tr_lvl                    = ', tr_lvl
     610          45 :          write(nu_diag,1010) ' tr_pond_cesm              = ', tr_pond_cesm
     611          45 :          write(nu_diag,1010) ' tr_pond_lvl               = ', tr_pond_lvl
     612          45 :          write(nu_diag,1010) ' tr_pond_topo              = ', tr_pond_topo
     613          45 :          write(nu_diag,1010) ' tr_aero                   = ', tr_aero
     614          45 :          write(nu_diag,1010) ' tr_fsd                    = ', tr_fsd
     615             : 
     616          45 :          nt_Tsfc = 1           ! index tracers, starting with Tsfc = 1
     617          45 :          ntrcr = 1             ! count tracers, starting with Tsfc = 1
     618             : 
     619          45 :          nt_qice = ntrcr + 1
     620          45 :          ntrcr = ntrcr + nilyr ! qice in nilyr layers
     621             : 
     622          45 :          nt_qsno = ntrcr + 1
     623          45 :          ntrcr = ntrcr + nslyr ! qsno in nslyr layers
     624             : 
     625          45 :          nt_sice = ntrcr + 1
     626          45 :          ntrcr = ntrcr + nilyr ! sice in nilyr layers
     627             : 
     628          45 :          nt_iage = max_ntrcr
     629          45 :          if (tr_iage) then
     630           2 :              ntrcr = ntrcr + 1
     631           2 :              nt_iage = ntrcr   ! chronological ice age
     632             :          endif
     633             : 
     634          45 :          nt_FY = max_ntrcr
     635          45 :          if (tr_FY) then
     636           2 :              ntrcr = ntrcr + 1
     637           2 :              nt_FY = ntrcr     ! area of first year ice
     638             :          endif
     639             : 
     640          45 :          nt_alvl = max_ntrcr
     641          45 :          nt_vlvl = max_ntrcr
     642          45 :          if (tr_lvl) then
     643          41 :              ntrcr = ntrcr + 1
     644          41 :              nt_alvl = ntrcr   ! area of level ice
     645          41 :              ntrcr = ntrcr + 1
     646          41 :              nt_vlvl = ntrcr   ! volume of level ice
     647             :          endif
     648             : 
     649          45 :          nt_apnd = max_ntrcr
     650          45 :          nt_hpnd = max_ntrcr
     651          45 :          nt_ipnd = max_ntrcr
     652          45 :          if (tr_pond) then            ! all explicit melt pond schemes
     653          36 :              ntrcr = ntrcr + 1
     654          36 :              nt_apnd = ntrcr
     655          36 :              ntrcr = ntrcr + 1
     656          36 :              nt_hpnd = ntrcr
     657          36 :              if (tr_pond_lvl) then
     658          29 :                  ntrcr = ntrcr + 1    ! refrozen pond ice lid thickness
     659          29 :                  nt_ipnd = ntrcr      ! on level-ice ponds (if frzpnd='hlid')
     660             :              endif
     661          36 :              if (tr_pond_topo) then
     662           5 :                  ntrcr = ntrcr + 1    ! 
     663           5 :                  nt_ipnd = ntrcr      ! refrozen pond ice lid thickness
     664             :              endif
     665             :          endif
     666             : 
     667          45 :          nt_fsd = max_ntrcr
     668          45 :          if (tr_fsd) then
     669           4 :              nt_fsd = ntrcr + 1       ! floe size distribution
     670           4 :              ntrcr = ntrcr + nfsd
     671             :          end if
     672             : 
     673          45 :          nt_isosno = max_ntrcr
     674          45 :          nt_isoice = max_ntrcr
     675          45 :          if (tr_iso) then             ! isotopes
     676           2 :              nt_isosno = ntrcr + 1
     677           2 :              ntrcr = ntrcr + n_iso    ! n_iso species in snow
     678           2 :              nt_isoice = ntrcr + 1
     679           2 :              ntrcr = ntrcr + n_iso    ! n_iso species in ice
     680             :          end if
     681             : 
     682          45 :          nt_aero = max_ntrcr - 4*n_aero
     683          45 :          if (tr_aero) then
     684           2 :              nt_aero = ntrcr + 1
     685           2 :              ntrcr = ntrcr + 4*n_aero ! 4 dEdd layers, n_aero species
     686             :          endif
     687             :               
     688          45 :          if (ntrcr > max_ntrcr-1) then
     689           0 :             write(nu_diag,*) 'max_ntrcr-1 < number of namelist tracers'
     690           0 :             write(nu_diag,*) 'max_ntrcr-1 = ',max_ntrcr-1,' ntrcr = ',ntrcr
     691           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     692             :          endif                               
     693             : 
     694          45 :          write(nu_diag,*) ' '
     695          45 :          write(nu_diag,1020) 'max_ntrcr = ', max_ntrcr
     696          45 :          write(nu_diag,1020) 'ntrcr = ', ntrcr
     697          45 :          write(nu_diag,*) ' '
     698          45 :          write(nu_diag,1020) 'nt_sice = ', nt_sice
     699          45 :          write(nu_diag,1020) 'nt_qice = ', nt_qice
     700          45 :          write(nu_diag,1020) 'nt_qsno = ', nt_qsno
     701          45 :          write(nu_diag,*)' '
     702          45 :          write(nu_diag,1020) 'ncat    = ', ncat
     703          45 :          write(nu_diag,1020) 'nilyr   = ', nilyr
     704          45 :          write(nu_diag,1020) 'nslyr   = ', nslyr
     705          45 :          write(nu_diag,1020) 'nblyr   = ', nblyr
     706          45 :          write(nu_diag,1020) 'nfsd    = ', nfsd
     707          45 :          write(nu_diag,1020) 'n_iso   = ', n_iso
     708          45 :          write(nu_diag,1020) 'n_aero  = ', n_aero
     709          45 :          write(nu_diag,*)' '
     710          45 :          write(nu_diag,1020) 'nx      = ', nx
     711          45 :          write(nu_diag,*)' '
     712             : 
     713             :  1000    format (a30,2x,f9.2)  ! a30 to align formatted, unformatted statements
     714             :  1005    format (a30,2x,f10.6)  ! float
     715             :  1010    format (a30,2x,l6)    ! logical
     716             :  1020    format (a30,2x,i6)    ! integer
     717             :  1030    format (a30,   a8)    ! character
     718             :  1040    format (a30,2x,6i6)   ! integer
     719             :  1050    format (a30,2x,6a6)   ! character
     720             : 
     721          45 :       if (formdrag) then
     722           3 :          if (nt_apnd==0) then
     723           0 :             write(nu_diag,*)'ERROR: nt_apnd:',nt_apnd
     724           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     725           3 :          elseif (nt_hpnd==0) then
     726           0 :             write(nu_diag,*)'ERROR: nt_hpnd:',nt_hpnd
     727           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     728           3 :          elseif (nt_ipnd==0) then
     729           0 :             write(nu_diag,*)'ERROR: nt_ipnd:',nt_ipnd
     730           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     731           3 :          elseif (nt_alvl==0) then
     732           0 :             write(nu_diag,*)'ERROR: nt_alvl:',nt_alvl
     733           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     734           3 :          elseif (nt_vlvl==0) then
     735           0 :             write(nu_diag,*)'ERROR: nt_vlvl:',nt_vlvl
     736           0 :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     737             :          endif
     738             :       endif
     739             : 
     740             :       !-----------------------------------------------------------------
     741             :       ! set Icepack values
     742             :       !-----------------------------------------------------------------
     743             : 
     744             :       call icepack_init_parameters(ustar_min_in=ustar_min, Cf_in=Cf, &
     745             :            albicev_in=albicev, albicei_in=albicei, ksno_in=ksno, &
     746             :            albsnowv_in=albsnowv, albsnowi_in=albsnowi, &
     747             :            natmiter_in=natmiter, ahmax_in=ahmax, shortwave_in=shortwave, &
     748             :            atmiter_conv_in = atmiter_conv, &
     749             :            albedo_type_in=albedo_type, R_ice_in=R_ice, R_pnd_in=R_pnd, &
     750             :            R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
     751             :            kstrength_in=kstrength, krdg_partic_in=krdg_partic, &
     752             :            krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &
     753             :            atmbndy_in=atmbndy, calc_strair_in=calc_strair, &
     754             :            formdrag_in=formdrag, highfreq_in=highfreq, &
     755             :            emissivity_in=emissivity, &
     756             :            kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, &
     757             :            dpscale_in=dpscale, frzpnd_in=frzpnd, &
     758             :            rfracmin_in=rfracmin, rfracmax_in=rfracmax, &
     759             :            pndaspect_in=pndaspect, hs1_in=hs1, hp1_in=hp1, &
     760             :            ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, &
     761             :            conduct_in=conduct, a_rapid_mode_in=a_rapid_mode, &
     762             :            Rac_rapid_mode_in=Rac_rapid_mode, &
     763             :            aspect_rapid_mode_in=aspect_rapid_mode, &
     764             :            dSdt_slow_mode_in=dSdt_slow_mode, &
     765             :            phi_c_slow_mode_in=phi_c_slow_mode, &
     766             :            phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
     767             :            tfrz_option_in=tfrz_option, kalg_in=kalg, &
     768             :            fbot_xfer_type_in=fbot_xfer_type, &
     769          45 :            wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec)
     770             :       call icepack_init_tracer_sizes(ntrcr_in=ntrcr, &
     771             :            ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, &
     772          45 :            nfsd_in=nfsd, n_iso_in=n_iso, n_aero_in=n_aero)
     773             :       call icepack_init_tracer_flags(tr_iage_in=tr_iage, &
     774             :            tr_FY_in=tr_FY, tr_lvl_in=tr_lvl, tr_aero_in=tr_aero, &
     775             :            tr_iso_in=tr_iso, &
     776             :            tr_pond_in=tr_pond, tr_pond_cesm_in=tr_pond_cesm, &
     777             :            tr_pond_lvl_in=tr_pond_lvl, &
     778          45 :            tr_pond_topo_in=tr_pond_topo, tr_fsd_in=tr_fsd)
     779             :       call icepack_init_tracer_indices(nt_Tsfc_in=nt_Tsfc, &
     780             :            nt_sice_in=nt_sice, nt_qice_in=nt_qice, &
     781             :            nt_qsno_in=nt_qsno, nt_iage_in=nt_iage, &
     782             :            nt_fy_in=nt_fy, nt_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl, &
     783             :            nt_apnd_in=nt_apnd, nt_hpnd_in=nt_hpnd, nt_ipnd_in=nt_ipnd, &
     784             :            nt_aero_in=nt_aero, nt_fsd_in=nt_fsd, &
     785          45 :            nt_isosno_in=nt_isosno, nt_isoice_in=nt_isoice)
     786             : 
     787          45 :       call icepack_warnings_flush(nu_diag)
     788          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     789           0 :           file=__FILE__,line= __LINE__)
     790             : 
     791          45 :       end subroutine input_data
     792             : 
     793             : !=======================================================================
     794             : 
     795             : ! Horizontal grid initialization:
     796             : !
     797             : ! author: Elizabeth C. Hunke, LANL
     798             : 
     799          45 :       subroutine init_grid2
     800             : 
     801             :       integer :: i
     802          16 :       real (kind=dbl_kind) :: pi, puny
     803             :       character(len=*), parameter :: subname='(init_grid2)'
     804             : 
     805             :       !-----------------------------------------------------------------
     806             :       ! query Icepack values
     807             :       !-----------------------------------------------------------------
     808             : 
     809          45 :       call icepack_query_parameters(pi_out=pi,puny_out=puny)
     810          45 :       call icepack_warnings_flush(nu_diag)
     811          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     812           0 :           file=__FILE__, line=__LINE__)
     813             : 
     814             :       !-----------------------------------------------------------------
     815             :       ! lat, lon, cell widths, angle, land mask
     816             :       !-----------------------------------------------------------------
     817             : 
     818         225 :       TLAT(:) = p5*pi  ! pi/2, North pole
     819          45 :       TLON(:) = c0
     820             : 
     821         180 :       do i = 2, nx
     822         180 :          TLAT(i) = TLAT(i-1) - p5*pi/180._dbl_kind ! half-deg increments
     823             :       enddo
     824             : 
     825         225 :       tmask(:) = .true.
     826          45 :       tmask(nx) = .false.  ! land in last grid cell
     827             : 
     828             :       !-----------------------------------------------------------------
     829             :       ! create hemisphere masks
     830             :       !-----------------------------------------------------------------
     831             : 
     832          45 :       lmask_n(:) = .false.
     833          45 :       lmask_s(:) = .false.
     834             : 
     835         225 :       do i = 1, nx
     836         180 :          if (TLAT(i) >= -puny) lmask_n(i) = .true. ! N. Hem.
     837         225 :          if (TLAT(i) <  -puny) lmask_s(i) = .true. ! S. Hem.
     838             :       enddo
     839             : 
     840          45 :       end subroutine init_grid2
     841             : 
     842             : !=======================================================================
     843             : 
     844             : ! Initialize state for the itd model
     845             : !
     846             : ! authors: C. M. Bitz, UW
     847             : !          William H. Lipscomb, LANL
     848             : 
     849          45 :       subroutine init_state
     850             : 
     851             :       use icepack_intfc, only: icepack_aggregate
     852             :       use icedrv_domain_size, only: ncat, nilyr, nslyr, nblyr, max_ntrcr
     853             :       use icedrv_domain_size, only: n_iso, n_aero, nfsd
     854             :       use icedrv_flux, only: sst, Tf, Tair, salinz, Tmltz
     855             :       use icedrv_state, only: trcr_depend, aicen, trcrn, vicen, vsnon
     856             :       use icedrv_state, only: aice0, aice, vice, vsno, trcr, aice_init
     857             :       use icedrv_state, only: n_trcr_strata, nt_strata, trcr_base
     858             : 
     859             :       integer (kind=int_kind) :: &
     860             :          i           , & ! horizontal indes
     861             :          k           , & ! vertical index
     862             :          it              ! tracer index
     863             : 
     864             :       logical (kind=log_kind) :: &
     865             :          heat_capacity   ! from icepack
     866             : 
     867             :       integer (kind=int_kind) :: ntrcr
     868             :       logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_fsd, tr_iso
     869             :       logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
     870             :       integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_fy
     871             :       integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, &
     872             :                                  nt_ipnd, nt_aero, nt_fsd, nt_isosno, nt_isoice
     873             : 
     874             :       character(len=*), parameter :: subname='(init_state)'
     875             : 
     876             :       !-----------------------------------------------------------------
     877             :       ! query Icepack values
     878             :       !-----------------------------------------------------------------
     879             : 
     880          45 :          call icepack_query_parameters(heat_capacity_out=heat_capacity)
     881          45 :          call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
     882             :          call icepack_query_tracer_flags(tr_iage_out=tr_iage, &
     883             :               tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, &
     884             :               tr_iso_out=tr_iso, &
     885             :               tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, &
     886          45 :               tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd)
     887             :          call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, &
     888             :               nt_sice_out=nt_sice, nt_qice_out=nt_qice, &
     889             :               nt_qsno_out=nt_qsno, nt_iage_out=nt_iage, nt_fy_out=nt_fy, &
     890             :               nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, &
     891             :               nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, &
     892             :               nt_ipnd_out=nt_ipnd, &
     893             :               nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
     894          45 :               nt_aero_out=nt_aero, nt_fsd_out=nt_fsd)
     895          45 :          call icepack_warnings_flush(nu_diag)
     896          45 :          if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     897           0 :              file=__FILE__,line= __LINE__)
     898             : 
     899             :       !-----------------------------------------------------------------
     900             :       ! Check number of layers in ice and snow.
     901             :       !-----------------------------------------------------------------
     902             :          if (nilyr < 1) then
     903             :             write (nu_diag,*) 'nilyr =', nilyr
     904             :             write (nu_diag,*) 'Must have at least one ice layer'
     905             :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     906             :          endif
     907             : 
     908             :          if (nslyr < 1) then
     909             :             write (nu_diag,*) 'nslyr =', nslyr
     910             :             write (nu_diag,*) 'Must have at least one snow layer'
     911             :             call icedrv_system_abort(file=__FILE__,line=__LINE__)
     912             :          endif
     913             : 
     914          45 :          if (.not.heat_capacity) then
     915             : 
     916           6 :             write (nu_diag,*) 'WARNING - Zero-layer thermodynamics'
     917             : 
     918             :             if (nilyr > 1) then
     919           0 :                write (nu_diag,*) 'nilyr =', nilyr
     920             :                write (nu_diag,*)        &
     921           0 :                     'Must have nilyr = 1 if ktherm = 0'
     922           0 :                call icedrv_system_abort(file=__FILE__,line=__LINE__)
     923             :             endif
     924             : 
     925             :             if (nslyr > 1) then
     926           0 :                write (nu_diag,*) 'nslyr =', nslyr
     927             :                write (nu_diag,*)        &
     928           0 :                     'Must have nslyr = 1 if heat_capacity = F'
     929           0 :                call icedrv_system_abort(file=__FILE__,line=__LINE__)
     930             :             endif
     931             : 
     932             :          endif   ! heat_capacity = F
     933             : 
     934             :       !-----------------------------------------------------------------
     935             :       ! Set tracer types
     936             :       !-----------------------------------------------------------------
     937             : 
     938          45 :       trcr_depend(nt_Tsfc) = 0 ! ice/snow surface temperature
     939         324 :       do k = 1, nilyr
     940         279 :          trcr_depend(nt_sice + k - 1) = 1 ! volume-weighted ice salinity
     941         324 :          trcr_depend(nt_qice + k - 1) = 1 ! volume-weighted ice enthalpy
     942             :       enddo
     943         102 :       do k = 1, nslyr
     944         102 :          trcr_depend(nt_qsno + k - 1) = 2 ! volume-weighted snow enthalpy
     945             :       enddo
     946          45 :       if (tr_iage) trcr_depend(nt_iage)  = 1   ! volume-weighted ice age
     947          45 :       if (tr_FY)   trcr_depend(nt_FY)    = 0   ! area-weighted first-year ice area
     948          45 :       if (tr_lvl)  trcr_depend(nt_alvl)  = 0   ! level ice area
     949          45 :       if (tr_lvl)  trcr_depend(nt_vlvl)  = 1   ! level ice volume
     950          45 :       if (tr_pond_cesm) then
     951           2 :                    trcr_depend(nt_apnd)  = 0           ! melt pond area
     952           2 :                    trcr_depend(nt_hpnd)  = 2+nt_apnd   ! melt pond depth
     953             :       endif
     954          45 :       if (tr_pond_lvl) then
     955          29 :                    trcr_depend(nt_apnd)  = 2+nt_alvl   ! melt pond area
     956          29 :                    trcr_depend(nt_hpnd)  = 2+nt_apnd   ! melt pond depth
     957          29 :                    trcr_depend(nt_ipnd)  = 2+nt_apnd   ! refrozen pond lid
     958             :       endif
     959          45 :       if (tr_pond_topo) then
     960           5 :                    trcr_depend(nt_apnd)  = 0           ! melt pond area
     961           5 :                    trcr_depend(nt_hpnd)  = 2+nt_apnd   ! melt pond depth
     962           5 :                    trcr_depend(nt_ipnd)  = 2+nt_apnd   ! refrozen pond lid
     963             :       endif
     964          45 :       if (tr_fsd) then
     965          41 :          do it = 1, nfsd
     966          41 :             trcr_depend(nt_fsd + it - 1) = 0    ! area-weighted floe size distribution
     967             :          enddo
     968             :       endif
     969          45 :       if (tr_iso) then  ! isotopes
     970           8 :          do it = 1, n_iso
     971           6 :             trcr_depend(nt_isosno) = 2          ! snow
     972           8 :             trcr_depend(nt_isoice) = 1          ! ice
     973             :          enddo
     974             :       endif
     975          45 :       if (tr_aero) then ! volume-weighted aerosols
     976           4 :          do it = 1, n_aero
     977           2 :             trcr_depend(nt_aero+(it-1)*4  ) = 2 ! snow
     978           2 :             trcr_depend(nt_aero+(it-1)*4+1) = 2 ! snow
     979           2 :             trcr_depend(nt_aero+(it-1)*4+2) = 1 ! ice
     980           4 :             trcr_depend(nt_aero+(it-1)*4+3) = 1 ! ice
     981             :          enddo
     982             :       endif
     983             : 
     984        1811 :       do it = 1, ntrcr
     985             :          ! mask for base quantity on which tracers are carried
     986        1766 :          if (trcr_depend(it) == 0) then      ! area
     987         156 :             trcr_base(it,1) = c1 
     988        1610 :          elseif (trcr_depend(it) == 1) then  ! ice volume
     989         612 :             trcr_base(it,2) = c1 
     990         998 :          elseif (trcr_depend(it) == 2) then  ! snow volume
     991         215 :             trcr_base(it,3) = c1 
     992             :          else
     993         783 :             trcr_base(it,1) = c1    ! default: ice area
     994         783 :             trcr_base(it,2) = c0
     995         783 :             trcr_base(it,3) = c0
     996             :          endif
     997             : 
     998             :          ! initialize number of underlying tracer layers
     999        1766 :          n_trcr_strata(it) = 0
    1000             :          ! default indices of underlying tracer layers
    1001        1766 :          nt_strata   (it,1) = 0
    1002        1811 :          nt_strata   (it,2) = 0
    1003             :       enddo
    1004             : 
    1005          45 :       if (tr_pond_cesm) then
    1006           2 :          n_trcr_strata(nt_hpnd)   = 1       ! melt pond depth
    1007           2 :          nt_strata    (nt_hpnd,1) = nt_apnd ! on melt pond area
    1008             :       endif
    1009          45 :       if (tr_pond_lvl) then
    1010          29 :          n_trcr_strata(nt_apnd)   = 1       ! melt pond area
    1011          29 :          nt_strata    (nt_apnd,1) = nt_alvl ! on level ice area
    1012          29 :          n_trcr_strata(nt_hpnd)   = 2       ! melt pond depth
    1013          29 :          nt_strata    (nt_hpnd,2) = nt_apnd ! on melt pond area
    1014          29 :          nt_strata    (nt_hpnd,1) = nt_alvl ! on level ice area
    1015          29 :          n_trcr_strata(nt_ipnd)   = 2       ! refrozen pond lid
    1016          29 :          nt_strata    (nt_ipnd,2) = nt_apnd ! on melt pond area
    1017          29 :          nt_strata    (nt_ipnd,1) = nt_alvl ! on level ice area
    1018             :       endif
    1019          45 :       if (tr_pond_topo) then
    1020           5 :          n_trcr_strata(nt_hpnd)   = 1       ! melt pond depth
    1021           5 :          nt_strata    (nt_hpnd,1) = nt_apnd ! on melt pond area
    1022           5 :          n_trcr_strata(nt_ipnd)   = 1       ! refrozen pond lid
    1023           5 :          nt_strata    (nt_ipnd,1) = nt_apnd ! on melt pond area
    1024             :       endif
    1025             : 
    1026             :       !-----------------------------------------------------------------
    1027             :       ! Set state variables
    1028             :       !-----------------------------------------------------------------
    1029             : 
    1030             :       call set_state_var (nx, &
    1031             :           Tair  (:),   sst  (:),     &
    1032             :           Tf    (:),                 &
    1033             :           salinz(:,:), Tmltz(:,:),   &
    1034             :           aicen (:,:), trcrn(:,:,:), &
    1035          45 :           vicen (:,:), vsnon(:,:))
    1036             : 
    1037             :       !-----------------------------------------------------------------
    1038             :       ! compute aggregate ice state and open water area
    1039             :       !-----------------------------------------------------------------
    1040             : 
    1041             :       !$OMP PARALLEL DO PRIVATE(i)
    1042         225 :       do i = 1, nx
    1043         180 :          aice(i) = c0
    1044         180 :          vice(i) = c0
    1045         180 :          vsno(i) = c0
    1046        9492 :          do it = 1, max_ntrcr
    1047        9492 :             trcr(i,it) = c0
    1048             :          enddo
    1049             : 
    1050         180 :          if (tmask(i)) &
    1051             :          call icepack_aggregate(ncat=ncat,                    &
    1052           0 :                                 trcrn=trcrn(i,1:ntrcr,:),     &
    1053             :                                 aicen=aicen(i,:),             &
    1054             :                                 vicen=vicen(i,:),             &
    1055             :                                 vsnon=vsnon(i,:),             &
    1056           0 :                                 trcr=trcr (i,1:ntrcr),        &
    1057          48 :                                 aice=aice (i),                &
    1058          48 :                                 vice=vice (i),                &
    1059          48 :                                 vsno=vsno (i),                &
    1060          48 :                                 aice0=aice0(i),               &
    1061             :                                 ntrcr=ntrcr,                  &
    1062           0 :                                 trcr_depend=trcr_depend(1:ntrcr),     &
    1063           0 :                                 trcr_base=trcr_base    (1:ntrcr,:),   &
    1064           0 :                                 n_trcr_strata=n_trcr_strata(1:ntrcr), &
    1065         231 :                                 nt_strata=nt_strata    (1:ntrcr,:))
    1066             : 
    1067         225 :          aice_init(i) = aice(i)
    1068             : 
    1069             :       enddo
    1070             :       !$OMP END PARALLEL DO
    1071             : 
    1072          45 :       call icepack_warnings_flush(nu_diag)
    1073          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1074           0 :           file=__FILE__, line=__LINE__)
    1075             : 
    1076          45 :       end subroutine init_state
    1077             : 
    1078             : !=======================================================================
    1079             : 
    1080             : ! Initialize state in each ice thickness category
    1081             : !
    1082             : ! authors: Elizabeth Hunke, LANL
    1083             : 
    1084          45 :       subroutine set_state_var (nx, &
    1085          90 :                                 Tair,     sst,    &
    1086          45 :                                 Tf,               &
    1087          90 :                                 salinz,   Tmltz,  &
    1088          90 :                                 aicen,    trcrn,  &
    1089          45 :                                 vicen,    vsnon)
    1090             : 
    1091             :       use icedrv_arrays_column, only: hin_max
    1092             :       use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, nfsd
    1093             :       use icedrv_arrays_column, only: floe_rad_c, floe_binwidth
    1094             : 
    1095             : 
    1096             :       integer (kind=int_kind), intent(in) :: &
    1097             :          nx          ! number of grid cells
    1098             : 
    1099             :       real (kind=dbl_kind), dimension (nx), intent(in) :: &
    1100             :          Tair       ! air temperature  (K)
    1101             : 
    1102             :       ! ocean values may be redefined here, unlike in CICE
    1103             :       real (kind=dbl_kind), dimension (nx), intent(inout) :: &
    1104             :          Tf     , & ! freezing temperature (C) 
    1105             :          sst        ! sea surface temperature (C) 
    1106             : 
    1107             :       real (kind=dbl_kind), dimension (nx,nilyr), &
    1108             :          intent(in) :: &
    1109             :          salinz , & ! initial salinity profile
    1110             :          Tmltz      ! initial melting temperature profile
    1111             : 
    1112             :       real (kind=dbl_kind), dimension (nx,ncat), &
    1113             :          intent(out) :: &
    1114             :          aicen , & ! concentration of ice
    1115             :          vicen , & ! volume per unit area of ice          (m)
    1116             :          vsnon     ! volume per unit area of snow         (m)
    1117             : 
    1118             :       real (kind=dbl_kind), dimension (nx,max_ntrcr,ncat), &
    1119             :          intent(out) :: &
    1120             :          trcrn     ! ice tracers
    1121             :                    ! 1: surface temperature of ice/snow (C)
    1122             : 
    1123             :       ! local variables
    1124             : 
    1125             :       integer (kind=int_kind) :: &
    1126             :          i     , & ! horizontal indices
    1127             :          k     , & ! ice layer index
    1128             :          n     , & ! thickness category index
    1129             :          it        ! tracer index
    1130             : 
    1131             :       real (kind=dbl_kind) :: &
    1132          32 :          Tsfc, sum, hbar, &
    1133          32 :          rhos, Lfresh, puny
    1134             : 
    1135             :       real (kind=dbl_kind), dimension(ncat) :: &
    1136         168 :          ainit, hinit    ! initial area, thickness
    1137             : 
    1138             :       real (kind=dbl_kind), dimension(nilyr) :: &
    1139         116 :          qin             ! ice enthalpy (J/m3)
    1140             : 
    1141             :       real (kind=dbl_kind), dimension(nslyr) :: &
    1142          36 :          qsn             ! snow enthalpy (J/m3)
    1143             : 
    1144             :       real (kind=dbl_kind), parameter :: &
    1145             :          hsno_init = 0.25_dbl_kind   ! initial snow thickness (m)
    1146             : 
    1147             :       logical (kind=log_kind) :: tr_brine, tr_lvl, tr_fsd
    1148             :       integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice, nt_fsd
    1149             :       integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl
    1150             : 
    1151             :       character(len=*), parameter :: subname='(set_state_var)'
    1152             : 
    1153             :       !-----------------------------------------------------------------
    1154             :       ! query Icepack values
    1155             :       !-----------------------------------------------------------------
    1156             : 
    1157             :       call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl, &
    1158          45 :         tr_fsd_out=tr_fsd)
    1159             :       call icepack_query_tracer_indices( nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
    1160             :            nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, &
    1161          45 :            nt_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl)
    1162          45 :       call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny)
    1163          45 :       call icepack_warnings_flush(nu_diag)
    1164          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1165           0 :          file=__FILE__,line= __LINE__)
    1166             : 
    1167             :       !-----------------------------------------------------------------
    1168             :       ! Initialize state variables.
    1169             :       ! If restarting, these values are overwritten.
    1170             :       !-----------------------------------------------------------------
    1171             : 
    1172         258 :       do n = 1, ncat
    1173        1065 :          do i = 1, nx
    1174         852 :             aicen(i,n) = c0
    1175         852 :             vicen(i,n) = c0
    1176         852 :             vsnon(i,n) = c0
    1177         852 :             trcrn(i,nt_Tsfc,n) = Tf(i)  ! surface temperature 
    1178             :             if (max_ntrcr >= 2) then
    1179       45744 :                do it = 2, max_ntrcr
    1180       45744 :                   trcrn(i,it,n) = c0
    1181             :                enddo
    1182             :             endif
    1183         852 :             if (tr_lvl)   trcrn(i,nt_alvl,n) = c1
    1184         852 :             if (tr_lvl)   trcrn(i,nt_vlvl,n) = c1
    1185         852 :             if (tr_brine) trcrn(i,nt_fbri,n) = c1
    1186        6384 :             do k = 1, nilyr
    1187        6384 :                trcrn(i,nt_sice+k-1,n) = salinz(i,k)
    1188             :             enddo
    1189        2157 :             do k = 1, nslyr
    1190        1944 :                trcrn(i,nt_qsno+k-1,n) = -rhos * Lfresh
    1191             :             enddo
    1192             :          enddo
    1193         213 :          ainit(n) = c0
    1194         258 :          hinit(n) = c0
    1195             :       enddo
    1196             : 
    1197             :       !-----------------------------------------------------------------
    1198             :       ! For Icepack testing, the grid vector is populated with several
    1199             :       ! different ice distributions, including ice-free, a single-
    1200             :       ! thickness slab, a full thickness distribution (as in CICE),
    1201             :       ! and land
    1202             :       !-----------------------------------------------------------------
    1203             : 
    1204          45 :       i = 1  ! ice-free
    1205             :              ! already initialized above 
    1206             : 
    1207             :       !-----------------------------------------------------------------
    1208             : 
    1209          45 :       i = 2  ! 2-m slab, no snow
    1210             :       if (3 <= ncat) then
    1211          42 :          n = 3
    1212          42 :          ainit(n) = c1  ! assumes we are using the default ITD boundaries
    1213          42 :          hinit(n) = c2
    1214             :       else
    1215           3 :          ainit(ncat) = c1
    1216           3 :          hinit(ncat) = c2
    1217             :       endif
    1218         258 :       do n = 1, ncat
    1219             :          ! ice volume, snow volume
    1220         213 :          aicen(i,n) = ainit(n)
    1221         213 :          vicen(i,n) = hinit(n) * ainit(n) ! m
    1222         213 :          vsnon(i,n) = c0
    1223             :          ! tracers
    1224          76 :          call icepack_init_trcr(Tair     = Tair(i),     &
    1225          76 :                                 Tf       = Tf(i),       &
    1226             :                                 Sprofile = salinz(i,:), &
    1227             :                                 Tprofile = Tmltz(i,:),  &
    1228             :                                 Tsfc     = Tsfc,        &
    1229             :                                 nilyr=nilyr, nslyr=nslyr, &
    1230         213 :                                 qin=qin(:), qsn=qsn(:))
    1231             : 
    1232             :          ! floe size distribution
    1233         213 :          if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
    1234             :                                   floe_rad_c=floe_rad_c,             &
    1235             :                                   floe_binwidth=floe_binwidth,       &
    1236          20 :                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
    1237             :          ! surface temperature
    1238         213 :          trcrn(i,nt_Tsfc,n) = Tsfc ! deg C
    1239             :          ! ice enthalpy, salinity 
    1240        1596 :          do k = 1, nilyr
    1241        1383 :             trcrn(i,nt_qice+k-1,n) = qin(k)
    1242        1596 :             trcrn(i,nt_sice+k-1,n) = salinz(i,k)
    1243             :          enddo
    1244             :          ! snow enthalpy
    1245         486 :          do k = 1, nslyr
    1246         486 :             trcrn(i,nt_qsno+k-1,n) = qsn(k)
    1247             :          enddo               ! nslyr
    1248             :          ! brine fraction
    1249         258 :          if (tr_brine) trcrn(i,nt_fbri,n) = c1
    1250             :       enddo                  ! ncat
    1251          45 :       call icepack_warnings_flush(nu_diag)
    1252          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1253           0 :           file=__FILE__, line=__LINE__)
    1254             :       
    1255             :       !-----------------------------------------------------------------
    1256             : 
    1257          45 :       i = 3  ! full thickness distribution
    1258             :       ! initial category areas in cells with ice
    1259          45 :       hbar = c3  ! initial ice thickness with greatest area
    1260             :       ! Note: the resulting average ice thickness 
    1261             :       ! tends to be less than hbar due to the
    1262             :       ! nonlinear distribution of ice thicknesses 
    1263             : 
    1264          45 :       sum = c0
    1265         258 :       do n = 1, ncat
    1266         213 :          if (n < ncat) then
    1267         168 :             hinit(n) = p5*(hin_max(n-1) + hin_max(n)) ! m
    1268             :          else                ! n=ncat
    1269          45 :             hinit(n) = (hin_max(n-1) + c1) ! m
    1270             :          endif
    1271             :          ! parabola, max at h=hbar, zero at h=0, 2*hbar
    1272         213 :          ainit(n) = max(c0, (c2*hbar*hinit(n) - hinit(n)**2))
    1273         258 :          sum = sum + ainit(n)
    1274             :       enddo
    1275         258 :       do n = 1, ncat
    1276         258 :          ainit(n) = ainit(n) / (sum + puny/ncat) ! normalize
    1277             :       enddo
    1278             :       
    1279         258 :       do n = 1, ncat
    1280             :          ! ice volume, snow volume
    1281         213 :          aicen(i,n) = ainit(n)
    1282         213 :          vicen(i,n) = hinit(n) * ainit(n) ! m
    1283         213 :          vsnon(i,n) = min(aicen(i,n)*hsno_init,p2*vicen(i,n))
    1284             :          ! tracers
    1285          76 :          call icepack_init_trcr(Tair     = Tair(i),     &
    1286          76 :                                 Tf       = Tf(i),       &
    1287             :                                 Sprofile = salinz(i,:), &
    1288             :                                 Tprofile = Tmltz(i,:),  &
    1289             :                                 Tsfc     = Tsfc,        &
    1290             :                                 nilyr=nilyr, nslyr=nslyr, &
    1291         213 :                                 qin=qin(:), qsn=qsn(:))
    1292             :          ! floe size distribution
    1293         213 :          if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
    1294             :                                   floe_rad_c=floe_rad_c,             &
    1295             :                                   floe_binwidth=floe_binwidth,       &
    1296          20 :                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
    1297             : 
    1298             :          ! surface temperature
    1299         213 :          trcrn(i,nt_Tsfc,n) = Tsfc ! deg C
    1300             :          ! ice enthalpy, salinity 
    1301        1596 :          do k = 1, nilyr
    1302        1383 :             trcrn(i,nt_qice+k-1,n) = qin(k)
    1303        1596 :             trcrn(i,nt_sice+k-1,n) = salinz(i,k)
    1304             :          enddo
    1305             :          ! snow enthalpy
    1306         486 :          do k = 1, nslyr
    1307         486 :             trcrn(i,nt_qsno+k-1,n) = qsn(k)
    1308             :          enddo               ! nslyr
    1309             :          ! brine fraction
    1310         258 :          if (tr_brine) trcrn(i,nt_fbri,n) = c1
    1311             :       enddo                  ! ncat
    1312          45 :       call icepack_warnings_flush(nu_diag)
    1313          45 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1314           0 :           file=__FILE__, line=__LINE__)
    1315             :       
    1316             :       !-----------------------------------------------------------------
    1317             :       
    1318             :       ! land
    1319             :       ! already initialized above (tmask = 0)
    1320          45 :       i = 4
    1321          45 :       sst(i) = c0
    1322          45 :       Tf(i) = c0
    1323             : 
    1324          45 :       end subroutine set_state_var
    1325             : 
    1326             : !=======================================================================
    1327             : 
    1328             : !  Initialize floe size distribution tracer (call prior to reading restart data)
    1329             : 
    1330          45 :       subroutine init_fsd
    1331             : 
    1332             :       use icedrv_arrays_column, only: wavefreq, dwavefreq, wave_sig_ht, &
    1333             :          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
    1334             :          d_afsd_wave, d_afsd_weld
    1335             : 
    1336          45 :       wavefreq       (:)   = c0
    1337          45 :       dwavefreq      (:)   = c0
    1338          45 :       wave_sig_ht    (:)   = c0
    1339          45 :       wave_spectrum  (:,:) = c0
    1340          45 :       d_afsd_newi    (:,:) = c0
    1341          45 :       d_afsd_latg    (:,:) = c0
    1342          45 :       d_afsd_latm    (:,:) = c0
    1343          45 :       d_afsd_wave    (:,:) = c0
    1344          45 :       d_afsd_weld    (:,:) = c0
    1345             : 
    1346          45 :       end subroutine init_fsd
    1347             : 
    1348             : !=======================================================================
    1349             : 
    1350             :   end module icedrv_init
    1351             : 
    1352             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd