MODULE diawri
  USE oce
  USE dom_oce
  USE phycst
  USE dianam
  USE diahth
  USE dynadv, ONLY: ln_dynadv_vec
  USE icb_oce
  USE icbdia
  USE ldftra
  USE ldfdyn
  USE sbc_oce
  USE sbc_ice
  USE sbcssr
  USE sbcwave
  USE wet_dry
  USE zdf_oce
  USE zdfdrg
  USE zdfmxl
  USE lbclnk
  USE in_out_manager
  USE diatmb
  USE dia25h
  USE iom
  USE ioipsl
  USE icewri
  USE lib_mpp
  USE timing
  USE diurnal_bulk
  USE cool_skin
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: dia_wri
  PUBLIC :: dia_wri_state
  PUBLIC :: dia_wri_alloc
  INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT
  INTEGER :: nb_T, ndim_bT
  INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU
  INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV
  INTEGER :: nid_W, nz_W, nh_W
  INTEGER :: ndex(1)
  INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
  INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
  INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
  CONTAINS
  INTEGER FUNCTION dia_wri_alloc()
    INTEGER, DIMENSION(2) :: ierr
    ierr = 0
    ALLOCATE(ndex_hT(jpi * jpj), ndex_T(jpi * jpj * jpk), ndex_hU(jpi * jpj), ndex_U(jpi * jpj * jpk), ndex_hV(jpi * jpj), ndex_V(jpi * jpj * jpk), STAT = ierr(1))
    dia_wri_alloc = MAXVAL(ierr)
    IF (lk_mpp) CALL mpp_sum(dia_wri_alloc)
  END FUNCTION dia_wri_alloc
  SUBROUTINE dia_wri(kt)
    INTEGER, INTENT( IN ) :: kt
    LOGICAL :: ll_print = .FALSE.
    CHARACTER(LEN = 40) :: clhstnam, clop, clmx
    INTEGER :: inum = 11
    INTEGER :: ji, jj, jk
    INTEGER :: ierr
    INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma
    INTEGER :: jn, ierror
    REAL(KIND = wp) :: zsto, zout, zmax, zjulian
    REAL(KIND = wp), DIMENSION(jpi, jpj) :: zw2d
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk) :: zw3d
    RETURN
    IF (ln_timing) CALL timing_start('dia_wri')
    IF (ninist == 1) THEN
      CALL dia_wri_state('output.init', kt)
      ninist = 0
    END IF
    ll_print = .FALSE.
    ll_print = ll_print .AND. lwp
    clop = "x"
    zsto = rdt
    clop = "ave(" // TRIM(clop) // ")"
    zout = nwrite * rdt
    zmax = (nitend - nit000 + 1) * rdt
    iimi = 1
    iima = jpi
    ijmi = 1
    ijma = jpj
    ipk = jpk
    it = kt
    itmod = kt - nit000 + 1
    IF (kt == nit000) THEN
      CALL ymds2ju(nyear, nmonth, nday, rdt, zjulian)
      zjulian = zjulian - adatrj
      IF (lwp) WRITE(numout, FMT = *)
      IF (lwp) WRITE(numout, FMT = *) 'Date 0 used :', nit000, ' YEAR ', nyear, ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
      IF (lwp) WRITE(numout, FMT = *) ' indexes of zoom = ', iimi, iima, ijmi, ijma, ' limit storage in depth = ', ipk
      IF (lwp) THEN
        CALL dia_nam(clhstnam, nwrite, ' ')
        CALL ctl_opn(inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', - 1, numout, lwp, narea)
        WRITE(inum, FMT = *) clhstnam
        CLOSE(UNIT = inum)
      END IF
      CALL dia_nam(clhstnam, nwrite, 'grid_T')
      IF (lwp) WRITE(numout, FMT = *) " Name of NETCDF file ", clhstnam
      CALL histbeg(clhstnam, jpi, glamt, jpj, gphit, iimi, iima - iimi + 1, ijmi, ijma - ijmi + 1, nit000 - 1, zjulian, rdt, nh_T, nid_T, domain_id = nidom, snc4chunks = snc4set)
      CALL histvert(nid_T, "deptht", "Vertical T levels", "m", ipk, gdept_1d, nz_T, "down")
      CALL wheneq(jpi * jpj * ipk, tmask, 1, 1., ndex_T, ndim_T)
      CALL wheneq(jpi * jpj, tmask, 1, 1., ndex_hT, ndim_hT)
      IF (ln_icebergs) THEN
        ALLOCATE(ndex_bT(jpi * jpj * nclasses), STAT = ierror)
        IF (lk_mpp) CALL mpp_sum(ierror)
        IF (ierror /= 0) THEN
          CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
          RETURN
        END IF
        CALL histvert(nid_T, "class", "Iceberg class", "number", nclasses, class_num, nb_T)
        ndim_bT = 3
        DO jn = 1, nclasses
          ndex_bT((jn - 1) * jpi * jpj + 1 : jn * jpi * jpj) = ndex_hT(1 : jpi * jpj)
        END DO
      END IF
      CALL dia_nam(clhstnam, nwrite, 'grid_U')
      IF (lwp) WRITE(numout, FMT = *) " Name of NETCDF file ", clhstnam
      CALL histbeg(clhstnam, jpi, glamu, jpj, gphiu, iimi, iima - iimi + 1, ijmi, ijma - ijmi + 1, nit000 - 1, zjulian, rdt, nh_U, nid_U, domain_id = nidom, snc4chunks = snc4set)
      CALL histvert(nid_U, "depthu", "Vertical U levels", "m", ipk, gdept_1d, nz_U, "down")
      CALL wheneq(jpi * jpj * ipk, umask, 1, 1., ndex_U, ndim_U)
      CALL wheneq(jpi * jpj, umask, 1, 1., ndex_hU, ndim_hU)
      CALL dia_nam(clhstnam, nwrite, 'grid_V')
      IF (lwp) WRITE(numout, FMT = *) " Name of NETCDF file ", clhstnam
      CALL histbeg(clhstnam, jpi, glamv, jpj, gphiv, iimi, iima - iimi + 1, ijmi, ijma - ijmi + 1, nit000 - 1, zjulian, rdt, nh_V, nid_V, domain_id = nidom, snc4chunks = snc4set)
      CALL histvert(nid_V, "depthv", "Vertical V levels", "m", ipk, gdept_1d, nz_V, "down")
      CALL wheneq(jpi * jpj * ipk, vmask, 1, 1., ndex_V, ndim_V)
      CALL wheneq(jpi * jpj, vmask, 1, 1., ndex_hV, ndim_hV)
      CALL dia_nam(clhstnam, nwrite, 'grid_W')
      IF (lwp) WRITE(numout, FMT = *) " Name of NETCDF file ", clhstnam
      CALL histbeg(clhstnam, jpi, glamt, jpj, gphit, iimi, iima - iimi + 1, ijmi, ijma - ijmi + 1, nit000 - 1, zjulian, rdt, nh_W, nid_W, domain_id = nidom, snc4chunks = snc4set)
      CALL histvert(nid_W, "depthw", "Vertical W levels", "m", ipk, gdepw_1d, nz_W, "down")
      CALL histdef(nid_T, "votemper", "Temperature", "C", jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout)
      CALL histdef(nid_T, "vosaline", "Salinity", "PSU", jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout)
      IF (.NOT. ln_linssh) THEN
        CALL histdef(nid_T, "vovvle3t", "Level thickness", "m", jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout)
        CALL histdef(nid_T, "vovvldep", "T point depth", "m", jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout)
        CALL histdef(nid_T, "vovvldef", "Squared level deformation", "%^2", jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout)
      END IF
      CALL histdef(nid_T, "sosstsst", "Sea Surface temperature", "C", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sosaline", "Sea Surface Salinity", "PSU", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sossheig", "Sea Surface Height", "m", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sowaflup", "Net Upward Water Flux", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sorunoff", "River runoffs", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sosfldow", "downward salt flux", "PSU/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      IF (ln_linssh) THEN
        CALL histdef(nid_T, "sosst_cd", "Concentration/Dilution term on temperature", "KgC/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "sosss_cd", "Concentration/Dilution term on salinity", "KgPSU/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      END IF
      CALL histdef(nid_T, "sohefldo", "Net Downward Heat Flux", "W/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "soshfldo", "Shortwave Radiation", "W/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "somixhgt", "Turbocline Depth", "m", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "somxl010", "Mixed Layer Depth 0.01", "m", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "soicecov", "Ice fraction", "[0,1]", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histdef(nid_T, "sowindsp", "wind speed at 10m", "m/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      IF (ln_icebergs) THEN
        CALL histdef(nid_T, "calving", "calving mass input", "kg/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "calving_heat", "calving heat flux", "XXXX", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "berg_floating_melt", "Melt rate of icebergs + bits", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "berg_stored_ice", "Accumulated ice mass by class", "kg", jpi, jpj, nh_T, nclasses, 1, nclasses, nb_T, 32, clop, zsto, zout)
        IF (ln_bergdia) THEN
          CALL histdef(nid_T, "berg_melt", "Melt rate of icebergs", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_buoy_melt", "Buoyancy component of iceberg melt rate", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_eros_melt", "Erosion component of iceberg melt rate", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_conv_melt", "Convective component of iceberg melt rate", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_virtual_area", "Virtual coverage by icebergs", "m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "bits_src", "Mass source of bergy bits", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "bits_melt", "Melt rate of bergy bits", "kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "bits_mass", "Bergy bit density field", "kg/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_mass", "Iceberg density field", "kg/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
          CALL histdef(nid_T, "berg_real_calving", "Calving into iceberg class", "kg/s", jpi, jpj, nh_T, nclasses, 1, nclasses, nb_T, 32, clop, zsto, zout)
        END IF
      END IF
      IF (.NOT. ln_cpl) THEN
        CALL histdef(nid_T, "sohefldp", "Surface Heat Flux: Damping", "W/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "sowafldp", "Surface Water Flux: Damping", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "sosafldp", "Surface salt flux: damping", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      END IF
      IF (ln_cpl .AND. nn_ice <= 1) THEN
        CALL histdef(nid_T, "sohefldp", "Surface Heat Flux: Damping", "W/m2", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "sowafldp", "Surface Water Flux: Damping", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
        CALL histdef(nid_T, "sosafldp", "Surface salt flux: Damping", "Kg/m2/s", jpi, jpj, nh_T, 1, 1, 1, - 99, 32, clop, zsto, zout)
      END IF
      clmx = "l_max(only(x))"
      CALL histend(nid_T, snc4chunks = snc4set)
      CALL histdef(nid_U, "vozocrtx", "Zonal Current", "m/s", jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout)
      IF (ln_wave .AND. ln_sdw) THEN
        CALL histdef(nid_U, "sdzocrtx", "Stokes Drift Zonal Current", "m/s", jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout)
      END IF
      CALL histdef(nid_U, "sozotaux", "Wind Stress along i-axis", "N/m2", jpi, jpj, nh_U, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histend(nid_U, snc4chunks = snc4set)
      CALL histdef(nid_V, "vomecrty", "Meridional Current", "m/s", jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout)
      IF (ln_wave .AND. ln_sdw) THEN
        CALL histdef(nid_V, "sdmecrty", "Stokes Drift Meridional Current", "m/s", jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout)
      END IF
      CALL histdef(nid_V, "sometauy", "Wind Stress along j-axis", "N/m2", jpi, jpj, nh_V, 1, 1, 1, - 99, 32, clop, zsto, zout)
      CALL histend(nid_V, snc4chunks = snc4set)
      CALL histdef(nid_W, "vovecrtz", "Vertical Velocity", "m/s", jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout)
      CALL histdef(nid_W, "votkeavt", "Vertical Eddy Diffusivity", "m2/s", jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout)
      CALL histdef(nid_W, "votkeavm", "Vertical Eddy Viscosity", "m2/s", jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout)
      IF (ln_zdfddm) THEN
        CALL histdef(nid_W, "voddmavs", "Salt Vertical Eddy Diffusivity", "m2/s", jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout)
      END IF
      IF (ln_wave .AND. ln_sdw) THEN
        CALL histdef(nid_W, "sdvecrtz", "Stokes Drift Vertical Current", "m/s", jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout)
      END IF
      CALL histend(nid_W, snc4chunks = snc4set)
      IF (lwp) WRITE(numout, FMT = *)
      IF (lwp) WRITE(numout, FMT = *) 'End of NetCDF Initialization'
      IF (ll_print) CALL FLUSH(numout)
    END IF
    IF (lwp .AND. MOD(itmod, nwrite) == 0) THEN
      WRITE(numout, FMT = *) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
      WRITE(numout, FMT = *) '~~~~~~ '
    END IF
    IF (.NOT. ln_linssh) THEN
      CALL histwrite(nid_T, "votemper", it, tsn(:, :, :, jp_tem) * e3t_n(:, :, :), ndim_T, ndex_T)
      CALL histwrite(nid_T, "vosaline", it, tsn(:, :, :, jp_sal) * e3t_n(:, :, :), ndim_T, ndex_T)
      CALL histwrite(nid_T, "sosstsst", it, tsn(:, :, 1, jp_tem) * e3t_n(:, :, 1), ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "sosaline", it, tsn(:, :, 1, jp_sal) * e3t_n(:, :, 1), ndim_hT, ndex_hT)
    ELSE
      CALL histwrite(nid_T, "votemper", it, tsn(:, :, :, jp_tem), ndim_T, ndex_T)
      CALL histwrite(nid_T, "vosaline", it, tsn(:, :, :, jp_sal), ndim_T, ndex_T)
      CALL histwrite(nid_T, "sosstsst", it, tsn(:, :, 1, jp_tem), ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "sosaline", it, tsn(:, :, 1, jp_sal), ndim_hT, ndex_hT)
    END IF
    IF (.NOT. ln_linssh) THEN
      !$ACC KERNELS
      zw3d(:, :, :) = ((e3t_n(:, :, :) - e3t_0(:, :, :)) / e3t_0(:, :, :) * 100 * tmask(:, :, :)) ** 2
      !$ACC END KERNELS
      CALL histwrite(nid_T, "vovvle3t", it, e3t_n(:, :, :), ndim_T, ndex_T)
      CALL histwrite(nid_T, "vovvldep", it, gdept_n(:, :, :), ndim_T, ndex_T)
      CALL histwrite(nid_T, "vovvldef", it, zw3d, ndim_T, ndex_T)
    END IF
    CALL histwrite(nid_T, "sossheig", it, sshn, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "sowaflup", it, (emp - rnf), ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "sorunoff", it, rnf, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "sosfldow", it, sfx, ndim_hT, ndex_hT)
    IF (ln_linssh) THEN
      !$ACC KERNELS
      zw2d(:, :) = emp(:, :) * tsn(:, :, 1, jp_tem)
      !$ACC END KERNELS
      CALL histwrite(nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT)
      !$ACC KERNELS
      zw2d(:, :) = emp(:, :) * tsn(:, :, 1, jp_sal)
      !$ACC END KERNELS
      CALL histwrite(nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT)
    END IF
    CALL histwrite(nid_T, "sohefldo", it, qns + qsr, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "soshfldo", it, qsr, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "somixhgt", it, hmld, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "somxl010", it, hmlp, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "soicecov", it, fr_i, ndim_hT, ndex_hT)
    CALL histwrite(nid_T, "sowindsp", it, wndm, ndim_hT, ndex_hT)
    IF (ln_icebergs) THEN
      CALL histwrite(nid_T, "calving", it, berg_grid % calving, ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "calving_heat", it, berg_grid % calving_hflx, ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "berg_floating_melt", it, berg_grid % floating_melt, ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "berg_stored_ice", it, berg_grid % stored_ice, ndim_bT, ndex_bT)
      IF (ln_bergdia) THEN
        CALL histwrite(nid_T, "berg_melt", it, berg_melt, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_buoy_melt", it, buoy_melt, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_eros_melt", it, eros_melt, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_conv_melt", it, conv_melt, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_virtual_area", it, virtual_area, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "bits_src", it, bits_src, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "bits_melt", it, bits_melt, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "bits_mass", it, bits_mass, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_mass", it, berg_mass, ndim_hT, ndex_hT)
        CALL histwrite(nid_T, "berg_real_calving", it, real_calving, ndim_bT, ndex_bT)
      END IF
    END IF
    IF (.NOT. ln_cpl) THEN
      CALL histwrite(nid_T, "sohefldp", it, qrp, ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "sowafldp", it, erp, ndim_hT, ndex_hT)
      IF (ln_ssr) zw2d(:, :) = erp(:, :) * tsn(:, :, 1, jp_sal) * tmask(:, :, 1)
      CALL histwrite(nid_T, "sosafldp", it, zw2d, ndim_hT, ndex_hT)
    END IF
    IF (ln_cpl .AND. nn_ice <= 1) THEN
      CALL histwrite(nid_T, "sohefldp", it, qrp, ndim_hT, ndex_hT)
      CALL histwrite(nid_T, "sowafldp", it, erp, ndim_hT, ndex_hT)
      IF (ln_ssr) zw2d(:, :) = erp(:, :) * tsn(:, :, 1, jp_sal) * tmask(:, :, 1)
      CALL histwrite(nid_T, "sosafldp", it, zw2d, ndim_hT, ndex_hT)
    END IF
    CALL histwrite(nid_U, "vozocrtx", it, un, ndim_U, ndex_U)
    CALL histwrite(nid_U, "sozotaux", it, utau, ndim_hU, ndex_hU)
    CALL histwrite(nid_V, "vomecrty", it, vn, ndim_V, ndex_V)
    CALL histwrite(nid_V, "sometauy", it, vtau, ndim_hV, ndex_hV)
    CALL histwrite(nid_W, "vovecrtz", it, wn, ndim_T, ndex_T)
    CALL histwrite(nid_W, "votkeavt", it, avt, ndim_T, ndex_T)
    CALL histwrite(nid_W, "votkeavm", it, avm, ndim_T, ndex_T)
    IF (ln_zdfddm) THEN
      CALL histwrite(nid_W, "voddmavs", it, avs, ndim_T, ndex_T)
    END IF
    IF (ln_wave .AND. ln_sdw) THEN
      CALL histwrite(nid_U, "sdzocrtx", it, usd, ndim_U, ndex_U)
      CALL histwrite(nid_V, "sdmecrty", it, vsd, ndim_V, ndex_V)
      CALL histwrite(nid_W, "sdvecrtz", it, wsd, ndim_T, ndex_T)
    END IF
    IF (kt == nitend) THEN
      CALL histclo(nid_T)
      CALL histclo(nid_U)
      CALL histclo(nid_V)
      CALL histclo(nid_W)
    END IF
    IF (ln_timing) CALL timing_stop('dia_wri')
  END SUBROUTINE dia_wri
  SUBROUTINE dia_wri_state(cdfile_name, kt)
    CHARACTER(LEN = *), INTENT( IN ) :: cdfile_name
    INTEGER, INTENT( IN ) :: kt
    CHARACTER(LEN = 32) :: clname
    CHARACTER(LEN = 40) :: clop
    INTEGER :: id_i, nz_i, nh_i
    INTEGER, DIMENSION(1) :: idex
    REAL(KIND = wp) :: zsto, zout, zmax, zjulian
    clname = cdfile_name
    IF (.NOT. Agrif_Root()) clname = TRIM(Agrif_CFixed()) // '_' // TRIM(clname)
    zsto = rdt
    clop = "inst(x)"
    zout = rdt
    zmax = (nitend - nit000 + 1) * rdt
    IF (lwp) WRITE(numout, FMT = *)
    IF (lwp) WRITE(numout, FMT = *) 'dia_wri_state : single instantaneous ocean state'
    IF (lwp) WRITE(numout, FMT = *) '~~~~~~~~~~~~~   and forcing fields file created '
    IF (lwp) WRITE(numout, FMT = *) '                and named :', clname, '.nc'
    CALL ymds2ju(nyear, nmonth, nday, rdt, zjulian)
    zjulian = zjulian - adatrj
    CALL histbeg(clname, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, nit000 - 1, zjulian, rdt, nh_i, id_i, domain_id = nidom, snc4chunks = snc4set)
    CALL histvert(id_i, "deptht", "Vertical T levels", "m", jpk, gdept_1d, nz_i, "down")
    CALL histdef(id_i, "vosaline", "Salinity", "PSU", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    CALL histdef(id_i, "votemper", "Temperature", "C", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    CALL histdef(id_i, "sossheig", "Sea Surface Height", "m", jpi, jpj, nh_i, 1, 1, 1, nz_i, 32, clop, zsto, zout)
    CALL histdef(id_i, "vozocrtx", "Zonal Current", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    CALL histdef(id_i, "vomecrty", "Meridional Current", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    CALL histdef(id_i, "vovecrtz", "Vertical Velocity", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    IF (ALLOCATED(ahtu)) THEN
      CALL histdef(id_i, "ahtu", "u-eddy diffusivity", "m2/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
      CALL histdef(id_i, "ahtv", "v-eddy diffusivity", "m2/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    END IF
    IF (ALLOCATED(ahmt)) THEN
      CALL histdef(id_i, "ahmt", "t-eddy viscosity", "m2/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
      CALL histdef(id_i, "ahmf", "f-eddy viscosity", "m2/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    END IF
    CALL histdef(id_i, "sowaflup", "Net Upward Water Flux", "Kg/m2/S", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    CALL histdef(id_i, "sohefldo", "Net Downward Heat Flux", "W/m2", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    CALL histdef(id_i, "soshfldo", "Shortwave Radiation", "W/m2", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    CALL histdef(id_i, "soicecov", "Ice fraction", "[0,1]", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    CALL histdef(id_i, "sozotaux", "Zonal Wind Stress", "N/m2", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    CALL histdef(id_i, "sometauy", "Meridional Wind Stress", "N/m2", jpi, jpj, nh_i, 1, 1, 1, - 99, 32, clop, zsto, zout)
    IF (.NOT. ln_linssh) THEN
      CALL histdef(id_i, "vovvldep", "T point depth", "m", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
      CALL histdef(id_i, "vovvle3t", "T point thickness", "m", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    END IF
    IF (ln_wave .AND. ln_sdw) THEN
      CALL histdef(id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
      CALL histdef(id_i, "sdmecrty", "Stokes Drift Merid", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
      CALL histdef(id_i, "sdvecrtz", "Stokes Drift Vert", "m/s", jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout)
    END IF
    IF (nn_ice == 2) THEN
      CALL ice_wri_state(kt, id_i, nh_i)
    END IF
    idex(1) = 1
    CALL histwrite(id_i, "votemper", kt, tsn(:, :, :, jp_tem), jpi * jpj * jpk, idex)
    CALL histwrite(id_i, "vosaline", kt, tsn(:, :, :, jp_sal), jpi * jpj * jpk, idex)
    CALL histwrite(id_i, "sossheig", kt, sshn, jpi * jpj, idex)
    CALL histwrite(id_i, "vozocrtx", kt, un, jpi * jpj * jpk, idex)
    CALL histwrite(id_i, "vomecrty", kt, vn, jpi * jpj * jpk, idex)
    CALL histwrite(id_i, "vovecrtz", kt, wn, jpi * jpj * jpk, idex)
    IF (ALLOCATED(ahtu)) THEN
      CALL histwrite(id_i, "ahtu", kt, ahtu, jpi * jpj * jpk, idex)
      CALL histwrite(id_i, "ahtv", kt, ahtv, jpi * jpj * jpk, idex)
    END IF
    IF (ALLOCATED(ahmt)) THEN
      CALL histwrite(id_i, "ahmt", kt, ahmt, jpi * jpj * jpk, idex)
      CALL histwrite(id_i, "ahmf", kt, ahmf, jpi * jpj * jpk, idex)
    END IF
    CALL histwrite(id_i, "sowaflup", kt, emp - rnf, jpi * jpj, idex)
    CALL histwrite(id_i, "sohefldo", kt, qsr + qns, jpi * jpj, idex)
    CALL histwrite(id_i, "soshfldo", kt, qsr, jpi * jpj, idex)
    CALL histwrite(id_i, "soicecov", kt, fr_i, jpi * jpj, idex)
    CALL histwrite(id_i, "sozotaux", kt, utau, jpi * jpj, idex)
    CALL histwrite(id_i, "sometauy", kt, vtau, jpi * jpj, idex)
    IF (.NOT. ln_linssh) THEN
      CALL histwrite(id_i, "vovvldep", kt, gdept_n(:, :, :), jpi * jpj * jpk, idex)
      CALL histwrite(id_i, "vovvle3t", kt, e3t_n(:, :, :), jpi * jpj * jpk, idex)
    END IF
    IF (ln_wave .AND. ln_sdw) THEN
      CALL histwrite(id_i, "sdzocrtx", kt, usd, jpi * jpj * jpk, idex)
      CALL histwrite(id_i, "sdmecrty", kt, vsd, jpi * jpj * jpk, idex)
      CALL histwrite(id_i, "sdvecrtz", kt, wsd, jpi * jpj * jpk, idex)
    END IF
    CALL histclo(id_i)
    IF (ninist /= 1) THEN
      CALL histclo(nid_T)
      CALL histclo(nid_U)
      CALL histclo(nid_V)
      CALL histclo(nid_W)
    END IF
  END SUBROUTINE dia_wri_state
END MODULE diawri


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Make damping matrices.
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0
!> @param[in] nn number of 1-D GLL nodes
!> @param[in] xq GLL nodes
!> @param[in] wq GLL weights
!> @param[in] dd spectral derivatives matrix
!> @param[in] rho mass density
!> @param[in] gamma damping coefficient
!> @param[in] A11 costant values for the bilinear map
!> @param[in] A12 costant values for the bilinear map
!> @param[in] A13 costant values for the bilinear map
!> @param[in] A21 costant values for the bilinear map
!> @param[in] A22 costant values for the bilinear map
!> @param[in] A23 costant values for the bilinear map
!> @param[in] A31 costant values for the bilinear map
!> @param[in] A32 costant values for the bilinear map
!> @param[in] A33 costant values for the bilinear map
!> @param[in] B11 costant values for the bilinear map
!> @param[in] B12 costant values for the bilinear map
!> @param[in] B13 costant values for the bilinear map
!> @param[in] B21 costant values for the bilinear map
!> @param[in] B22 costant values for the bilinear map
!> @param[in] B23 costant values for the bilinear map
!> @param[in] B31 costant values for the bilinear map
!> @param[in] B32 costant values for the bilinear map
!> @param[in] B33 costant values for the bilinear map
!> @param[in] GG1 costant values for the bilinear map
!> @param[in] GG2 costant values for the bilinear map
!> @param[in] GG3 costant values for the bilinear map
!> @param[in] DD1 costant values for the bilinear map
!> @param[in] DD2 costant values for the bilinear map
!> @param[in] DD3 costant values for the bilinear map
!> @param[out] mc_el viscous forces proportional to the velocity field
!> @param[out] mck_el viscous forces proportional to the displacement field

      subroutine MAKE_DAMPING_MATRIX(nn,ct,ww,dd,rho,gamma,&
                           A11,A12,A13,A21,A22,A23,&
                           A31,A32,A33,B11,B12,B13,&
                           B21,B22,B23,B31,B32,B33,&
                           GG1,GG2,GG3,DD1,DD2,DD3,&
                           mc_el,mck_el)

      implicit none

      integer*4 :: nn
      integer*4 :: i,j,k

      real*8 :: A11,A12,A13,A21,A22,A23,A31,A32,A33
      real*8 :: B11,B12,B13,B21,B22,B23,B31,B32,B33
      real*8 :: GG1,GG2,GG3,DD1,DD2,DD3
      real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
      real*8, dimension(nn) :: ct,ww

      real*8, dimension(nn,nn) :: dd

      real*8, dimension(nn,nn,nn) :: rho,gamma
      real*8, dimension(nn,nn,nn) :: mc_el,mck_el

      do k = 1,nn
         do j = 1,nn
            do i = 1,nn
               dxdx = A11 +B12*ct(k) +B13*ct(j) &
                    + GG1*ct(j)*ct(k)
               dydx = A21 +B22*ct(k) +B23*ct(j) &
                    + GG2*ct(j)*ct(k)
               dzdx = A31 +B32*ct(k) +B33*ct(j) &
                    + GG3*ct(j)*ct(k)

               dxdy = A12 +B11*ct(k) +B13*ct(i) &
                    + GG1*ct(k)*ct(i)
               dydy = A22 +B21*ct(k) +B23*ct(i) &
                    + GG2*ct(k)*ct(i)
               dzdy = A32 +B31*ct(k) +B33*ct(i) &
                    + GG3*ct(k)*ct(i)

               dxdz = A13 +B11*ct(j) +B12*ct(i) &
                    + GG1*ct(i)*ct(j)
               dydz = A23 +B21*ct(j) +B22*ct(i) &
                    + GG2*ct(i)*ct(j)
               dzdz = A33 +B31*ct(j) +B32*ct(i) &
                    + GG3*ct(i)*ct(j)

               det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
                     - dydz * (dxdx*dzdy - dzdx*dxdy) &
                     + dzdz * (dxdx*dydy - dydx*dxdy)

                    mc_el(i,j,k) = 2* gamma(i,j,k) * rho(i,j,k) &
                                 * det_j * ww(i) * ww(j) * ww(k)

                   mck_el(i,j,k) = (gamma(i,j,k)**2) * rho(i,j,k) &
                                * det_j * ww(i) * ww(j) * ww(k)

            enddo
         enddo
      enddo

      return

      end subroutine MAKE_DAMPING_MATRIX


! { dg-do run }
! pr18778 abort on endfile without opening unit
      program test
      implicit none
      integer i
      endfile(8)
      rewind(8)
      read(8,end=0023)i
      STOP 1! should never get here
      stop
 0023 continue
      close(8,status='delete')
      end


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE testNonLinearSolverInterface
USE IntrType
USE VectorTypes
USE MatrixTypes
USE NonLinearSolverModule

IMPLICIT NONE
PRIVATE

PUBLIC :: NonLinearSolverInterface_Test

!Sets up a test interface to find the roots of the following system of equations:
!  2*x1 + x1*x2 = 2
!  2*x2 - x1*x2**2 = 2
!The root occurs at x1 = 0.5, x2 = 2
TYPE,EXTENDS(NonLinearSolverInterface_Base) :: NonLinearSolverInterface_Test
  CONTAINS
    PROCEDURE,PASS :: eval => eval_test
ENDTYPE NonLinearSolverInterface_Test

CONTAINS
SUBROUTINE eval_test(this,x,y)
  CLASS(NonLinearSolverInterface_Test),INTENT(IN) :: this
  CLASS(VectorType),INTENT(INOUT) :: x
  CLASS(VectorType),INTENT(INOUT) :: y
  !
  REAL(SRK) :: xreal(2)

  CALL x%getAll(xreal)
  CALL y%set((/2.0_SRK*xreal(1) + xreal(1)*xreal(2) - 2.0_SRK, & ! 2*x1 + x1*x2 -2
      2.0_SRK*xreal(2) - xreal(1)*xreal(2)**2.0_SRK - 2.0_SRK/)) ! 2*x2 - x1*x2**2 -2

ENDSUBROUTINE eval_test
ENDMODULE testNonLinearSolverInterface
!
PROGRAM testNonLinearSolver
#include "UnitTest.h"
USE UnitTest
USE IntrType
USE ParameterLists
USE FutilityComputingEnvironmentModule
USE VectorTypes
USE MatrixTypes
USE LinearSolverTypes
USE NonLinearSolverModule
USE testNonLinearSolverInterface

IMPLICIT NONE

TYPE(FutilityComputingEnvironment),TARGET :: ce
TYPE(NonLinearSolverInterface_Test),TARGET :: ftest
TYPE(NonLinearSolver_Native) :: nativeSolver

CALL setupTest()
CREATE_TEST("testNonLinearSolver")

REGISTER_SUBTEST('Native, Newton''s Method',testNativeNewton)

FINALIZE_TEST()
CALL clearTest()
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE testNativeNewton()
  REAL(SRK) :: refsol(2)=(/0.5_SRK,2.0_SRK/),testsol(2)
  TYPE(ParamType) :: plist
  TYPE(RealVectorType) :: x
  TYPE(LinearSolverType_Direct) :: refLS
  CLASS(LinearSolverType_Base),ALLOCATABLE :: testLS

  CALL plist%add('NonLinearSolver -> n',2)
  CALL plist%add('NonLinearSolver -> method',NLSOLVER_METHOD_NEWTON)
  CALL plist%add('NonLinearSolver -> tolerance',1.0E-14_SRK)
  CALL nativeSolver%init(ce,ftest,plist)
  CALL plist%clear()

  COMPONENT_TEST('InitBase')
  ASSERT_EQ(testSolverMethod(nativeSolver),NLSOLVER_METHOD_NEWTON,'%solverMethod')
  ASSERT_EQ(testN(nativeSolver),2,'%n')
  ASSERT_EQ(testIterations(nativeSolver),-1,'%iterations')
  ASSERT_EQ(testTol(nativeSolver),1.0E-14_SRK,'%tol')
  ASSERT(ASSOCIATED(testCE(nativeSolver)),'%ce')
  ASSERT(ASSOCIATED(testFunc(nativeSolver)),'%func')

  COMPONENT_TEST('Init')
  ASSERT_EQ(testTPLType(nativeSolver),NLSOLVER_TPL_NATIVE,'%TPLType')
  CALL testLinSys(nativeSolver,testLS)
  ASSERT(testLS%isInit,'%linSys')
  ASSERT(SAME_TYPE_AS(refLS,testLS),'%linSys dynamic type')
  DEALLOCATE(testLS)
  ASSERT(nativeSolver%isInit,'%isInit')

  CALL plist%add('VectorType -> n',2)
  CALL x%init(plist)
  CALL plist%clear()

  COMPONENT_TEST('Solve')
  CALL x%set((/0.0_SRK,0.0_SRK/))
  CALL nativeSolver%solve(x)
  CALL x%get(testsol)
  ASSERT_APPROXEQ(testsol(1),refsol(1),'solution(1)')
  ASSERT_APPROXEQ(testsol(2),refsol(2),'solution(2)')
  ASSERT_EQ(testIterations(nativeSolver),7,'%iterations')
  CALL x%clear()

  CALL nativeSolver%clear()

  COMPONENT_TEST('ClearBase')
  ASSERT_EQ(testSolverMethod(nativeSolver),-1,'%solverMethod')
  ASSERT_EQ(testTPLType(nativeSolver),-1,'%TPLType')
  ASSERT_EQ(testN(nativeSolver),-1,'%n')
  ASSERT_EQ(testIterations(nativeSolver),-1,'%iterations')
  ASSERT_EQ(testTol(nativeSolver),1.0E-5_SRK,'%tol')
  CALL testLinSys(nativeSolver,testLS)
  ASSERT(.NOT.ALLOCATED(testLS),'%linSys')
  ASSERT(.NOT.ASSOCIATED(testCE(nativeSolver)),'%ce')
  ASSERT(.NOT.ASSOCIATED(testFunc(nativeSolver)),'ASSOCIATED %func')
  ASSERT(.NOT.nativeSolver%isInit,'%isInit')

ENDSUBROUTINE testNativeNewton
!
!-------------------------------------------------------------------------------
SUBROUTINE setupTest()

  ALLOCATE(ce%exceptHandler)
  ALLOCATE(ce%parEnv)

ENDSUBROUTINE setupTest
!
!-------------------------------------------------------------------------------
SUBROUTINE clearTest()

  DEALLOCATE(ce%exceptHandler)
  CALL ce%parEnv%world%finalize()
  DEALLOCATE(ce%parEnv)

ENDSUBROUTINE clearTest
!
ENDPROGRAM testNonLinearSolver

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
!      Driver for sounding ingest
       program snd_main

!      Steve Albers      May-1999       Original Version

       use mem_grid, ONLY: topo

!      Access global parameters
       character*8 c8_project
       character*9 a9_time

       call get_grid_dim_xy(ni,nj,istatus)
       if (istatus .ne. 1) then
           write (6,*) 'Error getting horizontal domain dimensions'
           stop
       endif

       call get_c8_project(c8_project,istatus)
       if (istatus .ne. 1) then
          write (6,*) 'Error getting c8_project'
          go to 999
       endif

       call get_laps_cycle_time(laps_cycle_time,istatus)
       if(istatus .ne. 1)stop

       call GETENV('LAPS_A9TIME',a9_time)
       call s_len(a9_time,ilen)

       if(ilen .eq. 9)then
           write(6,*)' systime (from env) = ',a9_time
           call i4time_fname_lp(a9_time,i4time_sys,istatus)
       else
           call get_systime(i4time_sys,a9_time,istatus)
           if(istatus .ne. 1)go to 999
           write(6,*)' systime = ',a9_time
       endif

       allocate( topo(ni,nj), STAT=istat_alloc )
       if(istat_alloc .ne. 0)then
           write(6,*)' ERROR: Could not allocate topo'
           stop
       endif

       call read_static_grid(ni,nj,'AVG',topo,istatus)
       if(istatus .ne. 1)then
           stop
       endif

       call snd_sub(ni,nj,c8_project,laps_cycle_time,i4time_sys)          ! I

 999   continue

       write(6,*)' Successful end of ingest_sounding program...'

       deallocate(topo)

       end

       subroutine snd_sub(ni,nj,c8_project,laps_cycle_time,i4time_sys)    ! I

       character*150 path_to_raw_raob,path_to_local_raob
       character*150 path_to_raw_drpsnd
       character*150 path_to_raw_satsnd
       character*150 path_to_raw_poessnd
       character*200 path_to_raw_tower
       character*200 path_to_raw_radiometer
       character*8 c8_raob_format, c8_drpsnd_format, c8_project
       character*13 filename13,cvt_i4time_wfo_fname13

       logical l_fill_ht

       l_fill_ht = .true. ! default value

       call get_snd_parms(path_to_raw_raob,path_to_local_raob
     1                   ,path_to_raw_drpsnd
     1                   ,path_to_raw_satsnd
     1                   ,path_to_raw_poessnd
     1                   ,path_to_raw_tower
     1                   ,path_to_raw_radiometer
     1                   ,l_fill_ht
     1                   ,istatus)

       if(istatus .ne. 1)goto999

       lun_out = 21

       if(c8_project .ne. 'RSA')then

!          Read RAOB data
           if(c8_project .eq. 'AFTAC')then
             c8_raob_format = 'NIMBUS'
           else
             c8_raob_format = c8_project
           endif
           write(6,*)
           write(6,*)' Call ingest_raob, format=',c8_raob_format
           call ingest_raob(path_to_raw_raob,c8_raob_format,i4time_sys
     1                     ,l_fill_ht,lun_out)

           call s_len(c8_project,lenc8)

!          Read Dropsonde data
           if(c8_project .eq. 'NIMBUS')then
               c8_drpsnd_format = c8_project

           elseif(c8_project .eq. 'CWB')then
               c8_drpsnd_format = c8_project

           elseif(c8_project(1:min(3,lenc8)) .eq. 'AIR')then
               c8_drpsnd_format = 'SND'

           else
               c8_drpsnd_format = 'AVAPS'

           endif

           write(6,*)
           write(6,*)' Call ingest_drpsnd, format=',c8_drpsnd_format
           call ingest_drpsnd(path_to_raw_drpsnd,c8_drpsnd_format
     1                       ,l_fill_ht,lun_out)

       else ! RSA project
           c8_raob_format = 'WFO'
           write(6,*)
           write(6,*)' Call ingest_raob, format=',c8_raob_format
           call ingest_raob(path_to_raw_raob,c8_raob_format,i4time_sys
     1                     ,l_fill_ht,lun_out)

           c8_raob_format = 'RSA'
           write(6,*)
           write(6,*)' Call ingest_raob, format=',c8_raob_format
           call ingest_raob(path_to_local_raob,c8_raob_format,i4time_sys
     1                     ,l_fill_ht,lun_out)

       endif ! RSA project

       write(6,*)
       write(6,*)' Call tower_driver_sub'
       maxobs = 6000
       itime_before = 900
       itime_after = 900

       call get_max_stations(maxsta, istatus)
       if(istatus .ne. 1)stop

       call tower_driver_sub(ni,nj,lun_out
     1                      ,maxobs,laps_cycle_time
     1                      ,path_to_raw_tower
     1                      ,itime_before,itime_after
     1                      ,maxsta
     1                      ,istatus)

!      Satellite soundings -- added goes (7/26/04 db)

       write(6,*)
       if(c8_project .eq. 'AFWA')then
           write(6,*)' Call ingest_satsnd for AFWA...'
           call ingest_satsnd(path_to_raw_satsnd,c8_project
     1                       ,i4time_sys,lun_out)
       else
           write(6,*)' Call ingest_satsnd for POES...'
           call ingest_satsnd(path_to_raw_poessnd,c8_project
     1                       ,i4time_sys,lun_out)

           write(6,*)
           write(6,*)' Call ingest_goessnd...'
           call ingest_goessnd (path_to_raw_satsnd, i4time_sys,
     1                          lun_out, istatus)
       endif

!      Radiometer data
       write(6,*)
       i4time_rad = (i4time_sys/3600) * 3600
 900   filename13 = cvt_i4time_wfo_fname13(i4time_rad)
       call s_len(path_to_raw_radiometer,len_path)
       if(len_path .gt. 0)then
           write(6,*)' Call get_radiometer_data, input file at '
     1              ,filename13
           call get_radiometer_data(i4time_sys,laps_cycle_time,ni,nj
     1                             ,i4time_sys - 855/2 ! Radiomtr cyc time / 2
     1                             ,i4time_sys + 855/2 ! Radiomtr cyc time / 2
     1                             ,path_to_raw_radiometer(1:len_path)
     1                                           //'/'//filename13
     1                             ,lun_out,istatus)
           if(istatus .ne. 1 .and. i4time_sys - i4time_rad .lt. 3600
     1                                                           )then
               write(6,*)' Trying radiometer for an hour earlier'
               i4time_rad = i4time_rad - 3600
               go to 900
           endif
       endif ! len_path > 0

       close(lun_out) ! Output SND file

 999   return

       end

       subroutine get_snd_parms(path_to_raw_raob,path_to_local_raob
     1                         ,path_to_raw_drpsnd
     1                         ,path_to_raw_satsnd
     1                         ,path_to_raw_poessnd
     1                         ,path_to_raw_tower
     1                         ,path_to_raw_radiometer
     1                         ,l_fill_ht
     1                         ,istatus)

       character*150 path_to_raw_raob,path_to_local_raob
     1              ,path_to_raw_drpsnd
     1              ,path_to_raw_satsnd
     1              ,path_to_raw_poessnd
     1              ,path_to_raw_radiometer

       character*200 path_to_raw_tower

       logical l_fill_ht

       namelist /snd_nl/ path_to_raw_raob,path_to_local_raob
     1                  ,path_to_raw_satsnd
     1                  ,path_to_raw_poessnd,path_to_raw_tower
     1                  ,path_to_raw_drpsnd,path_to_raw_radiometer
     1                  ,l_fill_ht

       character*150 static_dir,filename

       call get_directory('nest7grid',static_dir,len_dir)

       filename = static_dir(1:len_dir)//'/snd.nl'

       open(1,file=filename,status='old',err=900)
       read(1,snd_nl,err=901)
       close(1)

       istatus = 1
       return

  900  print*,'error opening file ',filename
       istatus = 0
       return

  901  print*,'error reading snd_nl in ',filename
       istatus = 0
       return

       end

      subroutine get_nlevels_snd(pressure_mb,height_m,r_missing_data
     1                          ,nlevels_raw,nlevels_snd)

      real pressure_mb(nlevels_raw)
      real height_m(nlevels_raw)

      do i = 1,nlevels_raw
          if(abs(pressure_mb(i)) .le. 2000. .or.
     1       abs(height_m(i)) .le. 1000000.)then
              nlevels_snd = i
          endif
      enddo ! i

      return
      end

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
great_rsd(m1,hh1).
great_rsd(k1,c1).
great_rsd(z1,l1).
great_rsd(p1,m1).
great_rsd(hh1,x1).
great_rsd(aa1,dd1).
great_rsd(a1,kk1).
great_rsd(cc1,g1).
great_rsd(cc1,jj1).
great_rsd(f1,jj1).
great_rsd(p1,w1).
great_rsd(l1,jj1).
great_rsd(ff1,c1).
great_rsd(t1,kk1).
great_rsd(f1,c1).
great_rsd(p1,x1).
great_rsd(m1,ii1).
great_rsd(t1,x1).
great_rsd(e1,a1).
great_rsd(ee1,m1).
great_rsd(cc1,v1).
great_rsd(cc1,ii1).
great_rsd(ff1,z1).
great_rsd(m1,kk1).
great_rsd(w1,hh1).
great_rsd(kk1,dd1).
great_rsd(g1,h1).
great_rsd(t1,c1).
great_rsd(w1,q1).
great_rsd(q1,jj1).
great_rsd(u1,ee1).
great_rsd(t1,dd1).
great_rsd(m1,k1).
great_rsd(f1,aa1).
great_rsd(v1,ee1).
great_rsd(c1,ii1).
great_rsd(p1,kk1).
great_rsd(e1,g1).
great_rsd(f1,hh1).
great_rsd(g1,c1).
great_rsd(u1,g1).
great_rsd(q1,h1).
great_rsd(z1,x1).
great_rsd(ff1,g1).
great_rsd(a1,ff1).
great_rsd(cc1,ee1).
great_rsd(g1,ii1).
great_rsd(p1,ii1).
great_rsd(e1,w1).
great_rsd(e1,ii1).
great_rsd(aa1,c1).
great_rsd(cc1,t1).
great_rsd(aa1,q1).
great_rsd(cc1,z1).
great_rsd(q1,z1).
great_rsd(e1,h1).
great_rsd(e1,m1).
great_rsd(e1,q1).
great_rsd(p1,v1).
great_rsd(l1,dd1).
great_rsd(ii1,jj1).
great_rsd(e1,jj1).
great_rsd(f1,dd1).
great_rsd(q1,g1).
great_rsd(e1,hh1).
great_rsd(t1,q1).
great_rsd(q1,c1).
great_rsd(ee1,l1).
great_rsd(aa1,ee1).
great_rsd(e1,x1).
great_rsd(m1,z1).
great_rsd(c1,dd1).
great_rsd(f1,k1).
great_rsd(ff1,l1).
great_rsd(cc1,q1).
great_rsd(q1,dd1).
great_rsd(v1,kk1).
great_rsd(ff1,hh1).
great_rsd(w1,jj1).
great_rsd(c1,l1).
great_rsd(a1,ee1).
great_rsd(w1,c1).
great_rsd(a1,dd1).
great_rsd(u1,z1).
great_rsd(c1,z1).
great_rsd(e1,t1).
great_rsd(kk1,hh1).
great_rsd(ii1,h1).
great_rsd(a1,l1).
great_rsd(a1,k1).
great_rsd(t1,l1).
great_rsd(v1,z1).
great_rsd(p1,hh1).
great_rsd(ff1,jj1).
great_rsd(g1,z1).
great_rsd(f1,q1).
great_rsd(f1,p1).
great_rsd(u1,aa1).
great_rsd(aa1,m1).
great_rsd(kk1,jj1).
great_rsd(v1,k1).
great_rsd(f1,a1).
great_rsd(u1,h1).
great_rsd(p1,g1).
great_rsd(a1,aa1).
great_rsd(dd1,x1).
great_rsd(p1,q1).
great_rsd(p1,z1).
great_rsd(f1,g1).
great_rsd(v1,g1).
great_rsd(aa1,jj1).
great_rsd(f1,h1).
great_rsd(cc1,hh1).
great_rsd(jj1,dd1).
great_rsd(kk1,z1).
great_rsd(z1,dd1).
great_rsd(p1,t1).
great_rsd(w1,u1).
great_rsd(p1,ee1).
great_rsd(kk1,x1).
great_rsd(a1,hh1).
great_rsd(u1,c1).
great_rsd(w1,m1).
great_rsd(p1,jj1).
great_rsd(w1,x1).
great_rsd(hh1,l1).
great_rsd(t1,z1).
great_rsd(t1,ff1).
great_rsd(ee1,kk1).
great_rsd(f1,ii1).
great_rsd(f1,v1).
great_rsd(k1,g1).
great_rsd(aa1,x1).
great_rsd(g1,dd1).
great_rsd(v1,t1).
great_rsd(q1,m1).
great_rsd(f1,z1).
great_rsd(w1,ff1).
great_rsd(p1,aa1).
great_rsd(e1,aa1).
great_rsd(t1,h1).
great_rsd(w1,dd1).
great_rsd(v1,m1).
great_rsd(w1,g1).
great_rsd(c1,x1).
great_rsd(k1,l1).
great_rsd(ff1,dd1).
great_rsd(f1,ff1).
great_rsd(e1,cc1).
great_rsd(aa1,z1).
great_rsd(cc1,u1).
great_rsd(c1,jj1).
great_rsd(e1,k1).
great_rsd(p1,u1).
great_rsd(f1,u1).
great_rsd(a1,h1).
great_rsd(ee1,h1).
great_rsd(q1,ee1).
great_rsd(e1,ff1).
great_rsd(kk1,l1).
great_rsd(cc1,ff1).
great_rsd(u1,l1).
great_rsd(p1,l1).
great_rsd(aa1,ff1).
great_rsd(ii1,l1).
great_rsd(l1,x1).
great_rsd(q1,l1).
great_rsd(h1,l1).
great_rsd(t1,k1).
great_rsd(cc1,dd1).
great_rsd(v1,aa1).
great_rsd(u1,ii1).
great_rsd(g1,x1).
great_rsd(a1,w1).
great_rsd(aa1,ii1).
great_rsd(q1,ii1).
great_rsd(cc1,h1).
great_rsd(v1,jj1).
great_rsd(kk1,h1).
great_rsd(t1,g1).
great_rsd(aa1,h1).
great_rsd(ee1,x1).
great_rsd(a1,cc1).
great_rsd(k1,x1).
great_rsd(w1,aa1).
great_rsd(u1,q1).
great_rsd(ee1,jj1).
great_rsd(v1,h1).
great_rsd(k1,h1).
great_rsd(e1,v1).
great_rsd(v1,x1).
great_rsd(aa1,kk1).
great_rsd(m1,g1).
great_rsd(m1,jj1).
great_rsd(m1,l1).
great_rsd(m1,c1).
great_rsd(aa1,k1).
great_rsd(u1,jj1).
great_rsd(u1,hh1).
great_rsd(a1,z1).
great_rsd(w1,t1).
great_rsd(h1,x1).
great_rsd(hh1,dd1).
great_rsd(t1,jj1).
great_rsd(m1,h1).
great_rsd(t1,m1).
great_rsd(p1,h1).
great_rsd(a1,ii1).
great_rsd(p1,c1).
great_rsd(ff1,x1).
great_rsd(w1,kk1).
great_rsd(hh1,jj1).
great_rsd(cc1,kk1).
great_rsd(g1,t1).
great_rsd(q1,w1).
great_rsd(c1,g1).
great_rsd(dd1,w1).
great_rsd(ii1,c1).
great_rsd(k1,w1).
great_rsd(q1,cc1).
great_rsd(t1,a1).
great_rsd(k1,p1).
great_rsd(h1,a1).
great_rsd(k1,ff1).
great_rsd(z1,k1).
great_rsd(jj1,m1).
great_rsd(q1,f1).
great_rsd(cc1,e1).
great_rsd(c1,u1).
great_rsd(t1,cc1).
great_rsd(aa1,t1).
great_rsd(ee1,p1).
great_rsd(l1,c1).
great_rsd(jj1,a1).
great_rsd(m1,e1).
great_rsd(z1,q1).
great_rsd(q1,aa1).
great_rsd(dd1,h1).
great_rsd(h1,c1).
great_rsd(h1,aa1).
great_rsd(dd1,a1).
great_rsd(hh1,w1).
great_rsd(dd1,hh1).
great_rsd(z1,w1).
great_rsd(ee1,t1).
great_rsd(jj1,g1).
great_rsd(ee1,q1).
great_rsd(c1,k1).
great_rsd(hh1,p1).
great_rsd(c1,q1).
great_rsd(x1,c1).
great_rsd(hh1,k1).
great_rsd(w1,a1).
great_rsd(l1,cc1).
great_rsd(ii1,cc1).
great_rsd(x1,p1).
great_rsd(ii1,q1).
great_rsd(g1,k1).
great_rsd(x1,t1).
great_rsd(hh1,t1).
great_rsd(x1,ff1).
great_rsd(k1,aa1).
great_rsd(aa1,e1).
great_rsd(ii1,a1).
great_rsd(g1,a1).
great_rsd(ii1,k1).
great_rsd(hh1,u1).
great_rsd(aa1,cc1).
great_rsd(ii1,e1).
great_rsd(dd1,m1).
great_rsd(ff1,p1).
great_rsd(v1,cc1).
great_rsd(hh1,aa1).
great_rsd(x1,z1).
great_rsd(dd1,q1).
great_rsd(l1,p1).
great_rsd(h1,e1).
great_rsd(c1,ff1).
great_rsd(kk1,e1).
great_rsd(c1,p1).
great_rsd(jj1,cc1).
great_rsd(m1,cc1).
great_rsd(jj1,ff1).
great_rsd(jj1,c1).
great_rsd(jj1,z1).
great_rsd(m1,t1).
great_rsd(h1,u1).
great_rsd(ee1,e1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!!$PROGRAM CRCPACKINGTEST
!!$
!!$  INTEGER                       :: N,DIM
!!$  INTEGER                       :: INFO
!!$  CHARACTER                     :: UPLO
!!$  COMPLEX*16,DIMENSION(5,5)     :: A
!!$
!!$  COMPLEX*16, DIMENSION(13) :: VALUES
!!$  INTEGER,    DIMENSION(13) :: COLUMNS
!!$  INTEGER,    DIMENSION(6) :: ROWINDEX
!!$
!!$  N    = SIZE(A,1)
!!$  DIM  =  SIZE(VALUES)
!!$  UPLO = "L"
!!$  A = 0
!!$  A(1,:) = (/ 1,-1, 0,-3, 0/)
!!$  A(2,:) = (/-2, 5, 0, 0, 0/)
!!$  A(3,:) = (/ 0, 0, 4, 6, 4/)
!!$  A(4,:) = (/-4, 0, 2, 0, 0/)
!!$  A(5,:) = (/ 0, 8, 0, 0,-5/)
!!$
!!$  VALUES   = 0
!!$  COLUMNS  = 0
!!$  ROWINDEX = 0
!!$  INFO     = 0
!!$
!!$  CALL VARCRCPACKING(N,DIM,UPLO,A,VALUES,COLUMNS,ROWINDEX,INFO)
!!$  CALL WRITE_MATRIX(N,N, REAL(A))
!!$  !WRITE(*,*) VALUES
!!$  CALL WRITE_MATRIX(1,DIM,REAL(VALUES))
!!$  CALL WRITE_MATRIX_INT(1,DIM,COLUMNS)
!!$  CALL WRITE_MATRIX_INT(1,N+1,ROWINDEX)
!!$
!!$END PROGRAM CRCPACKINGTEST

!!$SUBROUTINE WRITE_MATRIX(N1,N2,A)
!!$  DOUBLE PRECISION, DIMENSION(N1,N2) :: A
!!$  CHARACTER(LEN=55) STRING
!!$  CHARACTER(LEN=55) aux_char
!!$  integer :: aux
!!$
!!$  aux = int(UBOUND(A,2))
!!$  write(aux_char,"(I4)") aux
!!$  aux_char = trim(aux_char)
!!$  write(string,"(A1,I4,A6)") "(",aux,"E15.6)"
!!$  !!string = trim(string)
!!$  !WRITE(*,*) string, aux,LBOUND(A,1),UBOUND(A,1),LBOUND(A,2),UBOUND(A,2)
!!$
!!$  DO I = LBOUND(A,1), UBOUND(A,1)
!!$     WRITE(*,string) (A(I,J), J = LBOUND(A,2), UBOUND(A,2))
!!$  END DO
!!$  WRITE(*,*)
!!$!  WRITE(*,*)
!!$  !101 format(string)
!!$END SUBROUTINE WRITE_MATRIX
!!$
!!$
!!$SUBROUTINE WRITE_MATRIX_INT(N1,N2,A)
!!$  INTEGER, DIMENSION(N1,N2) :: A
!!$  WRITE(*,*)
!!$  DO I = LBOUND(A,1), UBOUND(A,1)
!!$     WRITE(*,*) (A(I,J), J = LBOUND(A,2), UBOUND(A,2))
!!$  END DO
!!$END SUBROUTINE WRITE_MATRIX_INT
!!$

SUBROUTINE VARCRCPACKING(N,DIM,UPLO,zero,A,VALUES,COLUMNS,ROWINDEX,INFO)

  IMPLICIT NONE
  INTEGER,                   INTENT(IN)    :: N
  INTEGER,                   INTENT(INOUT) :: INFO,DIM
  CHARACTER,                 INTENT(IN)    :: UPLO
  DOUBLE PRECISION,          INTENT(IN)    :: ZERO
  COMPLEX*16,DIMENSION(N,N), INTENT(IN)    :: A

  COMPLEX*16, DIMENSION(DIM), INTENT(OUT) :: VALUES
  INTEGER,    DIMENSION(DIM), INTENT(OUT) :: COLUMNS
  INTEGER,    DIMENSION(N+1), INTENT(OUT) :: ROWINDEX

  INTEGER I,J, counter,aux

  VALUES   = 0
  counter  =  1
  columns  = -1
  ROWINDEX = -1

  SELECT CASE (UPLO)

  CASE('F')
     DO I=1,N
        DO J=1,N
           IF(ABS(A(I,J)).GT.ZERO .OR. I.EQ.J ) THEN
              VALUES(counter)  = A(I,J)
              COLUMNS(counter) = J
              IF(ROWINDEX(I).EQ.-1) ROWINDEX(I) = counter
              counter = counter + 1
           END IF
        END DO
        IF(I.EQ.N-1) aux = COUNTER
     END DO
     AUX = COUNTER - 1 - AUX
     ROWINDEX(N+1) = aux + ROWINDEX(N) + 1

  CASE ('U')
     DO I=1,N
        DO J=I,N
           IF(ABS(A(I,J)).GT.ZERO .OR. I.EQ.J) THEN
              VALUES(counter)  = A(I,J)
              COLUMNS(counter) = J
              IF(ROWINDEX(I).EQ.-1) ROWINDEX(I) = counter
              counter = counter + 1
           END IF
        END DO
        IF(I.EQ.N-1) aux =COUNTER
     END DO
     AUX = COUNTER - 1 - AUX

     ROWINDEX(N+1) = aux + ROWINDEX(N) + 1

  CASE('L')
     DO I=1,N
        DO J=1,I
           IF(ABS(A(I,J)).GT.ZERO  .OR. I.EQ.J) THEN
              VALUES(counter)  = A(I,J)
              COLUMNS(counter) = J
              IF(ROWINDEX(I).EQ.-1) ROWINDEX(I) = counter
              counter = counter + 1
           END IF
        END DO
        IF(I.EQ.N-1) aux = COUNTER
     END DO
     AUX = COUNTER - 1 - AUX

     ROWINDEX(N+1) = aux + ROWINDEX(N) + 1

  END SELECT

  DIM = counter - 1

END SUBROUTINE VARCRCPACKING

max_stars_repo_name: gsinuco/OPENMMF

max_issues_repo_name: gsinuco/OPENMMF

max_forks_repo_name: gsinuco/OPENMMF


text:
      subroutine read_radar_v01_cdf(filename,nxv01,nyv01,nzv01,
     &cref_com,cvel_com,cnyq_com,refd,veld,nyqd,istatus)
c
c note: nxv01= lon
c       nyv01= lat
c       nzv01= level
c
      real*4      refd(nxv01,nyv01,nzv01)
      real*4      veld(nxv01,nyv01,nzv01)
      real*4      nyqd(nxv01,nyv01,nzv01)

      character       cvel_com(nzv01)*126
      character       cref_com(nzv01)*126
      character       cnyq_com(nzv01)*126
      character*(*)   filename

      include 'netcdf.inc'
C
C  Open netcdf File for reading
C
      n=index(filename,' ')-1
      print*,'Open netCDF file',filename(1:n)
      nf_status = NF_OPEN(filename,NF_NOWRITE,nf_fid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'NF_OPEN ',filename(1:n)
      endif
C
C  Fill all dimension values
C
C
C Get size of lat
C
      nf_status = NF_INQ_DIMID(nf_fid,'lat',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim lat'
      endif
      nf_status = NF_INQ_DIMLEN(nf_fid,nf_vid,lat)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim lat'
      endif
C
C Get size of level
C
      nf_status = NF_INQ_DIMID(nf_fid,'level',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim level'
      endif
      nf_status = NF_INQ_DIMLEN(nf_fid,nf_vid,level)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim level'
      endif
C
C Get size of lon
C
      nf_status = NF_INQ_DIMID(nf_fid,'lon',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim lon'
      endif
      nf_status = NF_INQ_DIMLEN(nf_fid,nf_vid,lon)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim lon'
      endif
C
C Get size of n_2d_grids
C
      nf_status = NF_INQ_DIMID(nf_fid,'n_2d_grids',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim n_2d_grids'
      endif
      nf_status = NF_INQ_DIMLEN(nf_fid,nf_vid,n_2d_grids)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'dim n_2d_grids'
      endif
      if(lon.ne.nxv01)then
         print *,'lon ne nxv01'
         stop
      endif
      if(lat.ne.nyv01)then
         print *,'lat ne nyv01'
         stop
      endif
      if(level.ne.nzv01)then
         print *,'level ne nzv01'
         stop
      endif

      call main_sub(nf_fid , nyv01, nzv01, nxv01, n_2d_grids,
     +cref_com,cvel_com,cnyq_com,refd,veld,nyqd)

      end
C
C
      subroutine main_sub(nf_fid, lat, level, lon, n_2d_grids,
     +refd_comment,veld_comment,nyqd_comment,refd,veld,nyqd)
      include 'netcdf.inc'
      integer lat, level, lon, n_2d_grids, nf_fid, nf_vid, nf_status
      character*18 asctime
      character*12 laps_domain_file
      character*132 model
      character*126 nyqd_comment(level)
      character*132 origin
      character*126 refd_comment(level)
      character*126 veld_comment(level)
      integer fctimes_var, imax, jmax, kdim, kmax, level_var(level),
     +     lvl(n_2d_grids), num_variables, nyqd_fcinv(level),
     +     refd_fcinv(level), veld_fcinv(level), version

      real nyqd( lon,  lat, level), refd( lon,  lat, level), veld(
     +     lon,  lat, level)

      call read_netcdf(nf_fid , lat, level, lon, n_2d_grids, asctime,
     +     fctimes_var, imax, jmax, kdim, kmax, laps_domain_file,
     +     level_var, lvl, model, num_variables, nyqd, nyqd_comment,
     +     nyqd_fcinv, origin, refd, refd_comment, refd_fcinv, veld,
     +     veld_comment, veld_fcinv, version)

C
C The netcdf variables are filled - your code goes here
C
      return
      end
      subroutine read_netcdf(nf_fid , lat, level, lon, n_2d_grids,
     +     asctime, fctimes_var, imax, jmax, kdim, kmax,
     +     laps_domain_file, level_var, lvl, model, num_variables,
     +     nyqd, nyqd_comment, nyqd_fcinv, origin, refd,
     +     refd_comment, refd_fcinv, veld, veld_comment, veld_fcinv,
     +     version)

      include 'netcdf.inc'
      integer lat, level, lon, n_2d_grids, nf_fid, nf_vid, nf_status

      character*18 asctime
      character*12 laps_domain_file
      character*132 model
      character*126 nyqd_comment(level)
      character*132 origin
      character*126 refd_comment(level)
      character*126 veld_comment(level)
      integer fctimes_var, imax, jmax, kdim, kmax, level_var(level),
     +     lvl(n_2d_grids), num_variables, nyqd_fcinv(level),
     +     refd_fcinv(level), veld_fcinv(level), version

      real nyqd( lon,  lat, level), refd( lon,  lat, level), veld(
     +     lon,  lat, level)

C
C     Variable        NETCDF Long Name
C      asctime
        nf_status = NF_INQ_VARID(nf_fid,'asctime',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var asctime'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,asctime)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ asctime '
      endif
C
C     Variable        NETCDF Long Name
C      laps_domain_file
        nf_status = NF_INQ_VARID(nf_fid,'laps_domain_file',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var laps_domain_file'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,laps_domain_file)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ laps_domain_file '
      endif
C
C     Variable        NETCDF Long Name
C      model
        nf_status = NF_INQ_VARID(nf_fid,'model',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var model'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,model)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ model '
      endif
C
C     Variable        NETCDF Long Name
C      nyqd_comment
       nf_status = NF_INQ_VARID(nf_fid,'nyqd_comment',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var nyqd_comment'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,nyqd_comment)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ nyqd_comment '
      endif
C
C     Variable        NETCDF Long Name
C      origin
        nf_status = NF_INQ_VARID(nf_fid,'origin',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var origin'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,origin)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ origin '
      endif
C
C     Variable        NETCDF Long Name
C      refd_comment
        nf_status = NF_INQ_VARID(nf_fid,'refd_comment',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var refd_comment'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,refd_comment)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ refd_comment '
      endif
C
C     Variable        NETCDF Long Name
C      veld_comment
        nf_status = NF_INQ_VARID(nf_fid,'veld_comment',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var veld_comment'
      endif
        nf_status = NF_GET_VAR_TEXT(nf_fid,nf_vid,veld_comment)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ veld_comment '
      endif
C
C     Variable        NETCDF Long Name
C      fctimes_var  "forecast times"
C
        nf_status = NF_INQ_VARID(nf_fid,'fctimes',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var fctimes'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,fctimes_var)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ fctimes_var '
      endif
C
C     Variable        NETCDF Long Name
C      imax
        nf_status = NF_INQ_VARID(nf_fid,'imax',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var imax'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,imax)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ imax '
      endif
C
C     Variable        NETCDF Long Name
C      jmax
        nf_status = NF_INQ_VARID(nf_fid,'jmax',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var jmax'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,jmax)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ jmax '
      endif
C
C     Variable        NETCDF Long Name
C      kdim
        nf_status = NF_INQ_VARID(nf_fid,'kdim',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var kdim'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,kdim)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ kdim '
      endif
C
C     Variable        NETCDF Long Name
C      kmax
        nf_status = NF_INQ_VARID(nf_fid,'kmax',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var kmax'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,kmax)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ kmax '
      endif
C
C     Variable        NETCDF Long Name
C      level_var    "level"
C
        nf_status = NF_INQ_VARID(nf_fid,'level',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var level'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,level_var)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ level_var '
      endif
C
C     Variable        NETCDF Long Name
C      lvl
        nf_status = NF_INQ_VARID(nf_fid,'lvl',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var lvl'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,lvl)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ lvl '
      endif
C
C     Variable        NETCDF Long Name
C      num_variables
        nf_status = NF_INQ_VARID(nf_fid,'num_variables',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var num_variables'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,num_variables)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ num_variables '
      endif
C
C     Variable        NETCDF Long Name
C      nyqd_fcinv
        nf_status = NF_INQ_VARID(nf_fid,'nyqd_fcinv',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var nyqd_fcinv'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,nyqd_fcinv)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ nyqd_fcinv '
      endif
C
C     Variable        NETCDF Long Name
C      refd_fcinv
        nf_status = NF_INQ_VARID(nf_fid,'refd_fcinv',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var refd_fcinv'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,refd_fcinv)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ refd_fcinv '
      endif
C
C     Variable        NETCDF Long Name
C      veld_fcinv
        nf_status = NF_INQ_VARID(nf_fid,'veld_fcinv',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var veld_fcinv'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,veld_fcinv)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ veld_fcinv '
      endif
C
C     Variable        NETCDF Long Name
C      version
        nf_status = NF_INQ_VARID(nf_fid,'version',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var version'
      endif
        nf_status = NF_GET_VAR_INT(nf_fid,nf_vid,version)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ version '
      endif
C
C     Variable        NETCDF Long Name
C      nyqd         "3D radar nyquist velocity"
C
      print*,'Read nyqd'
      nf_status = NF_INQ_VARID(nf_fid,'nyqd',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var nyqd'
      endif
        nf_status = NF_GET_VAR_REAL(nf_fid,nf_vid,nyqd)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ nyqd '
      endif
C
C     Variable        NETCDF Long Name
C      refd         "3D radar"
C
      print*,'Read refd'
      nf_status = NF_INQ_VARID(nf_fid,'refd',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var refd'
      endif
        nf_status = NF_GET_VAR_REAL(nf_fid,nf_vid,refd)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ refd '
      endif
C
C     Variable        NETCDF Long Name
C      veld         "3D radar"
C
      print*,'Read veld'
      nf_status = NF_INQ_VARID(nf_fid,'veld',nf_vid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in var veld'
      endif
        nf_status = NF_GET_VAR_REAL(nf_fid,nf_vid,veld)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'in NF_GET_VAR_ veld '
      endif
      nf_status = nf_close(nf_fid)
      if(nf_status.ne.NF_NOERR) then
        print *, NF_STRERROR(nf_status)
        print *,'nf_close'
      endif

      return
      end


text:
!*****************************************************************************************
!> author: Jacob Williams
!  date: 12/22/2015
!  license: BSD
!
!  Root-finding routines. Used by [[ddeabm_module]].

    module root_module

    use ddeabm_kinds

    implicit none

    private

    !real parameters:
    real(wp),parameter :: zero  = 0.0_wp
    real(wp),parameter :: half  = 0.5_wp
    real(wp),parameter :: one   = 1.0_wp
    real(wp),parameter :: two   = 2.0_wp
    real(wp),parameter :: three = 3.0_wp
    real(wp),parameter :: eps   = epsilon(one)  !! original code had d1mach(4)

    abstract interface
        function func(t) result(f)  !! interface to the [[zeroin]] input function
        import :: wp
        implicit none
        real(wp),intent(in)  :: t  !! Independant variable for the function.
        real(wp)             :: f  !! The function evaluated at `t`.
        end function func
    end interface

    public :: zeroin      ! main routine

    contains
!*****************************************************************************************

!*****************************************************************************************
!>
!  Find a zero of the function \( f(x) \) in the given interval
!  \( [a_x,b_x] \) to within a tolerance \( 4 \epsilon |x| + tol \),
!  where \( \epsilon \) is the relative machine precision defined as
!  the smallest representable number such that \( 1.0 + \epsilon > 1.0 \).
!
!  It is assumed that \( f(a_x) \) and \( f(b_x) \) have opposite signs.
!
!# References
!  * R. P. Brent, "[An algorithm with guaranteed convergence for
!    finding a zero of a function](http://maths-people.anu.edu.au/~brent/pd/rpb005.pdf)",
!    The Computer Journal, Vol 14, No. 4., 1971.
!  * R. P. Brent, "[Algorithms for minimization without derivatives](http://maths-people.anu.edu.au/~brent/pub/pub011.html)",
!    Prentice-Hall, Inc., 1973.
!
!# See also
!  1. [zeroin.f](http://www.netlib.org/go/zeroin.f) from Netlib

    subroutine zeroin(f,ax,bx,tol,xzero,fzero,iflag,fax,fbx)

    use iso_fortran_env, only: error_unit

    implicit none

    procedure(func)              :: f       !! the function to find the root of
    real(wp),intent(in)          :: ax      !! left endpoint of initial interval
    real(wp),intent(in)          :: bx      !! right endpoint of initial interval
    real(wp),intent(in)          :: tol     !! desired length of the interval of uncertainty of the final result (>=0)
    real(wp),intent(out)         :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out)         :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)          :: iflag   !! status flag (`-1`=error, `0`=root found)
    real(wp),intent(in),optional :: fax     !! if `f(ax)` is already known, it can be input here
    real(wp),intent(in),optional :: fbx     !! if `f(ax)` is already known, it can be input here

    real(wp) :: a,b,c,d,e,fa,fb,fc,tol1,xm,p,q,r,s

    tol1 = eps+one

    a=ax
    b=bx

    if (present(fax)) then
        fa = fax
    else
        fa=f(a)
    end if
    if (present(fbx)) then
        fb = fbx
    else
        fb=f(b)
    end if

    !check trivial cases first:
    if (fa==zero) then

        iflag = 0
        xzero = a
        fzero = fa

    else if (fb==zero) then

        iflag = 0
        xzero = b
        fzero = fb

    else if (fa*(fb/abs(fb))<zero) then  ! check that f(ax) and f(bx) have different signs

        c=a
        fc=fa
        d=b-a
        e=d

        do

            if (abs(fc)<abs(fb)) then
                a=b
                b=c
                c=a
                fa=fb
                fb=fc
                fc=fa
            end if

            tol1 = two*eps*abs(b)+half*tol
            xm = half*(c-b)
            if ((abs(xm)<=tol1).or.(fb==zero)) exit

            ! see if a bisection is forced
            if ((abs(e)>=tol1).and.(abs(fa)>abs(fb))) then
                s=fb/fa
                if (a/=c) then
                    ! inverse quadratic interpolation
                    q=fa/fc
                    r=fb/fc
                    p=s*(two*xm*q*(q-r)-(b-a)*(r-one))
                    q=(q-one)*(r-one)*(s-one)
                else
                    ! linear interpolation
                    p=two*xm*s
                    q=one-s
                end if
                if (p<=zero) then
                    p=-p
                else
                    q=-q
                end if
                s=e
                e=d
                if (((two*p)>=(three*xm*q-abs(tol1*q))) .or. &
                    (p>=abs(half*s*q))) then
                    d=xm
                    e=d
                else
                    d=p/q
                end if
            else
                d=xm
                e=d
            end if

            a=b
            fa=fb
            if (abs(d)<=tol1) then
                if (xm<=zero) then
                    b=b-tol1
                else
                    b=b+tol1
                end if
            else
                b=b+d
            end if
            fb=f(b)
            if ((fb*(fc/abs(fc)))>zero) then
                c=a
                fc=fa
                d=b-a
                e=d
            end if

        end do

        iflag = 0
        xzero = b
        fzero = fb

    else

        iflag = -1
        write(error_unit,'(A)')&
            'Error in zeroin: f(ax) and f(bx) do not have different signs.'

    end if

    end subroutine zeroin
!*****************************************************************************************

!*****************************************************************************************
    end module root_module
!*****************************************************************************************

max_stars_repo_name: jacobwilliams/ddeabm


text:
      subroutine particle_setup(np,x_start,y_start,z_start,
     1     time_start,cell_start,x1,y1,z1,ibout_last,
     2     index,index_np_cell,adv_start)
c s kelkar March 28 06
c set up initial variables to start march along a particle trajectory
c********************************************************************
c INPUT VARIABLES
c     np : particle #
c
c OUTOUT VARIABLES
c     cell_start :
c     x_start, y_start,z_start,
c     x1,y1,z1
c     iin,jin,kin
c     time_start
c     ibout_last
c     nx,ny,nz
c     index_np_cell
c     cell_index_sg(np)=index_np_cell
c INTERNAL VARIABLES
c     cell_index(np)
c     cell_packed(cell_index(np))
c     nxyz(ip_touched_cells(cell_start),1:3):
c     x_packed, y_packed, z_packed (cell_index(np)) :
c     time_packed(cell_index(np)) :
c
c********************************************************************

      use comparttr
      use comparttr_sg

      implicit none

      integer np,cell_start,ibout_last
      integer index,index_np_cell, i_sg_in
      integer nx,ny,nz,iin,jin,kin

      real*8 x_start, y_start,z_start,time_start,adv_start
      real*8 x1,y1,z1,xin,yin,zin

      time_start=time_packed(cell_index(np))
      x_start=x_packed(cell_index(np))
      y_start=y_packed(cell_index(np))
      z_start=z_packed(cell_index(np))
      cell_start=cell_packed(cell_index(np))
      adv_start=adv_packed(cell_index(np))
      ibout_last=0
      call global_local(cell_start,x_start,y_start,z_start,
     1     x1,y1,z1)

      nx=nxyz(ip_touched_cells(cell_start),1)
      ny=nxyz(ip_touched_cells(cell_start),2)
      nz=nxyz(ip_touched_cells(cell_start),3)
      call xyz_ijk(cell_start,nx,ny,nz,x1,y1,z1,iin,jin,kin)
c      call local_global(id_cell,x,y,z,xg,yg,zg)
      call path_arrays_sg(np,cell_start,x_start, y_start,z_start,
     1     time_start,iin,jin,kin,i_sg_in,index,index_np_cell,
     2     adv_start)

      cell_index_sg(np)=index_np_cell

      return

      end

c..................................................................

max_stars_repo_name: lanl/PLUMECALC

max_issues_repo_name: lanl/PLUMECALC

max_forks_repo_name: lanl/PLUMECALC


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE SOLID_STRAIN_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OGEL, WRITE_FEMAP )

! Calculates element specific strain outputs from array STRESS (generated in subr ELEM_STRE_STRN_ARRAYS) for solid elements

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  SOLID_STRAIN_OUTPUTS_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO
      USE MODEL_STUF, ONLY            :  STRAIN, TYPE
      USE CC_OUTPUT_DESCRIBERS, ONLY  :  STRN_OPT
      USE LINK9_STUFF, ONLY           :  OGEL
      USE FEMAP_ARRAYS, ONLY          :  FEMAP_EL_VECS
      USE PARAMS, ONLY                :  POST

      USE SOLID_STRAIN_OUTPUTS_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'SOLID_STRAIN_OUTPUTS'
      CHARACTER( 1*BYTE), INTENT(IN)  :: WRITE_OGEL         ! If 'Y' then write data to array OGEL
      CHARACTER( 1*BYTE), INTENT(IN)  :: WRITE_FEMAP        ! If 'Y' then write data to array FEMAP_EL_VECS

      INTEGER(LONG), INTENT(IN)       :: SIZE_ALLOCATED     ! No. of rows allocated to array that will be written to
      INTEGER(LONG), INTENT(IN)       :: NUM_FEMAP_ROWS     ! Number of rows that will be written to FEMAP arrays
      INTEGER(LONG), INTENT(INOUT)    :: NUM1               ! Cum rows written to OGEL prior to running this subr
      INTEGER(LONG)                   :: J                  ! DO loop indices
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = SOLID_STRAIN_OUTPUTS_BEGEND

      REAL(DOUBLE)                    :: MEAN               ! Mean strains
      REAL(DOUBLE)                    :: PRINCIPAL_STRAIN(3)! Principal strains
      REAL(DOUBLE)                    :: SIG_OCT            ! Octrahedral normal strain
      REAL(DOUBLE)                    :: TAU_OCT            ! Octrahedral shear  strain
      REAL(DOUBLE)                    :: VONMISES           ! von Mises strain

      INTRINSIC DMAX1,DMIN1

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Calc engineering strains from array STRAIN and put into array OGEL

      IF ((TYPE(1:4) == 'HEXA') .OR. (TYPE(1:5) == 'PENTA') .OR. (TYPE(1:5) == 'TETRA')) THEN
         CALL PRINCIPAL_3D ( STRAIN, PRINCIPAL_STRAIN, MEAN, VONMISES, SIG_OCT, TAU_OCT )
         IF (WRITE_OGEL == 'Y') THEN
            NUM1 = NUM1 + 1
            IF (NUM1 > SIZE_ALLOCATED) THEN
               WRITE(ERR,9200) SUBR_NAME,SIZE_ALLOCATED
               WRITE(F06,9200) SUBR_NAME,SIZE_ALLOCATED
               FATAL_ERR = FATAL_ERR + 1
               CALL OUTA_HERE ( 'Y' )                      ! Coding error (dim of array OGEL too small), so quit
            ENDIF
            DO J=1,6
               OGEL(NUM1,J) = STRAIN(J)
            ENDDO
            IF (STRN_OPT == 'VONMISES') THEN
               OGEL(NUM1,7) = VONMISES
               OGEL(NUM1,8) = ZERO
            ELSE
               OGEL(NUM1,7) = SIG_OCT
               OGEL(NUM1,8) = TAU_OCT
            ENDIF
         ENDIF
         IF ((POST /= 0) .AND. (WRITE_FEMAP == 'Y')) THEN
            DO J=1,6
               FEMAP_EL_VECS(NUM_FEMAP_ROWS,J) = STRAIN(J)
            ENDDO
            FEMAP_EL_VECS(NUM_FEMAP_ROWS, 7) = PRINCIPAL_STRAIN(1)
            FEMAP_EL_VECS(NUM_FEMAP_ROWS, 8) = PRINCIPAL_STRAIN(2)
            FEMAP_EL_VECS(NUM_FEMAP_ROWS, 9) = PRINCIPAL_STRAIN(3)
            FEMAP_EL_VECS(NUM_FEMAP_ROWS,10) = MEAN
            IF (STRN_OPT == 'VONMISES') THEN
               FEMAP_EL_VECS(NUM_FEMAP_ROWS,11) = VONMISES
               FEMAP_EL_VECS(NUM_FEMAP_ROWS,12) = ZERO
            ELSE
               FEMAP_EL_VECS(NUM_FEMAP_ROWS,11) = SIG_OCT
               FEMAP_EL_VECS(NUM_FEMAP_ROWS,12) = TAU_OCT
            ENDIF
         ENDIF

      ELSE

         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,9203) SUBR_NAME, TYPE
         WRITE(F06,9203) SUBR_NAME, TYPE
         CALL OUTA_HERE ( 'Y' )

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 9200 FORMAT(' *ERROR  9200: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' ARRAY OGEL WAS ALLOCATED TO HAVE ',I12,' ROWS. ATTEMPT TO WRITE TO OGEL BEYOND THIS')

 9203 FORMAT(' *ERROR  9203: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' INCORRECT ELEMENT TYPE = "',A,'"')

! **********************************************************************************************************************************

      END SUBROUTINE SOLID_STRAIN_OUTPUTS

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


text: Bambu specializes in Vietnamese Chè, sweet beverages that are consumed as dessert.  They also have boba milk teas, hot/cold/blended coffee drinks, and smoothies.  From February 20February 22, theyre running a Grand Opening promotional price of $3 for all drinks, smoothies, and desserts.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine plane4(co,node,nodep,a,b,c,d)
!
!     given are 5 nodes: node and nodep(1) up to nodep(4)
!
!     first, the node among the nodep's with the largest
!     distance from node is eliminated. The remaining nodes
!     are stored into noden(1) up to noden(3).
!
!     Then, the equation of the plane through the
!     nodes noden(1),noden(2) and noden(3) in the form
!     a*x+b*y+c*z+d=0 such that the triangle through the
!     nodes noden(1),noden(2),nopen(3) is numbered clockwise
!     when looking in the direction of vector (a,b,c)
!
      implicit none
!
      integer nodep(4),i,j,noden(3),node,kflag,idist(4),n
!
      real*8 co(3,*),a,b,c,d,dd,p12(3),p23(3),p31(3),dist(4)
!
      kflag=2
      n=4
!
!     determining the distance of the nodep's to node
!
      do i=1,4
         dist(i)=((co(1,nodep(i))-co(1,node))**2+
     &            (co(2,nodep(i))-co(2,node))**2+
     &            (co(3,nodep(i))-co(3,node))**2)
         idist(i)=nodep(i)
c         write(*,*) nodep(i),dist(i)
      enddo
!
!     sorting the distances
!
      call dsort(dist,idist,n,kflag)
!
!     storing the 3 closest nodes in noden
!
      j=0
      do i=1,4
         if(nodep(i).eq.idist(4)) cycle
         j=j+1
         noden(j)=nodep(i)
      enddo
c      write(*,*) 'noden ',(noden(i),i=1,3)
!
!     sides of the triangle
!
      do i=1,3
         p12(i)=co(i,noden(2))-co(i,noden(1))
         p23(i)=co(i,noden(3))-co(i,noden(2))
         p31(i)=co(i,noden(1))-co(i,noden(3))
      enddo
!
!     normalized vector normal to the triangle: xn = p12 x p23
!
      a=p12(2)*p23(3)-p12(3)*p23(2)
      b=p12(3)*p23(1)-p12(1)*p23(3)
      c=p12(1)*p23(2)-p12(2)*p23(1)
      dd=dsqrt(a*a+b*b+c*c)
      a=a/dd
      b=b/dd
      c=c/dd
!
!     determining the inhomogeneous term
!
      d=-a*co(1,noden(1))-b*co(2,noden(1))-c*co(3,noden(1))
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop