! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + + ! + glide_diagnostics.f90 - part of the Glimmer-CISM ice model + ! + + ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010 ! Glimmer-CISM contributors - see AUTHORS file for list of contributors ! ! This file is part of Glimmer-CISM. ! ! Glimmer-CISM 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, either version 2 of the License, or (at ! your option) any later version. ! ! Glimmer-CISM 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 Glimmer-CISM. If not, see . ! ! Glimmer-CISM is hosted on BerliOS.de: ! https://developer.berlios.de/projects/glimmer-cism/ ! ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ module glide_diagnostics !*FD subroutines for computing various useful diagnostics ! Author: William Lipscomb, LANL contains subroutine glide_write_diag (model, time, idiag, jdiag) !*FD Writes diagnostic output use glimmer_log use glimmer_paramets, only: thk0, len0, vel0, tim0, GLC_DEBUG use glimmer_physcon, only: scyr use glide_types implicit none type(glide_global_type), intent(in) :: model ! model instance real(dp), intent(in) :: time ! current time in years integer, intent(in), optional :: idiag, jdiag real(dp) :: & tot_area, & ! total ice area (km^2) tot_volume, & ! total ice volume (km^3) tot_energy, & ! total ice energy (J) tot_age, & ! sum of volume*age max_thck, & ! max ice thickness (m) min_temp, & ! min ice temperature (deg C) mean_thck, & ! mean ice thickness (m) mean_temp, & ! mean ice temperature (deg C) mean_age, & ! mean ice age (yr) max_spd_sfc, & ! max surface ice speed (m/yr) max_spd_bas, & ! max basal ice speed (m/yr) thck, & ! thickness spd, & ! speed age ! age integer :: i, j, k, & imin, jmin, kmin, & imax, jmax, kmax, & ewn, nsn, upn ! model%numerics%ewn, etc. character(len=100) :: message real(dp), parameter :: & eps = 1.0e-11_dp ! small number ewn = model%general%ewn nsn = model%general%nsn upn = model%general%upn !----------------------------------------------------------------- ! Compute and write global diagnostics !----------------------------------------------------------------- print*, ' ' print*, 'Writing diagnostics to log file, time(yr) =', time call write_log(' ') write(message,'(a32,f24.16)') 'Global diagnostic output, time =', time call write_log(trim(message), type = GM_DIAGNOSTIC) call write_log(' ') ! total ice area tot_area = 0._dp do j = 1, nsn do i = 1, ewn if (model%geometry%thck(i,j) > model%numerics%thklim) then tot_area = tot_area + 1.0_dp endif enddo enddo tot_area = tot_area * model%numerics%dew * model%numerics%dns * len0**2 write(message,'(a25,e24.16)') 'Total ice area (km^2) ', & tot_area*1.0e-6_dp ! convert to km^2 call write_log(trim(message), type = GM_DIAGNOSTIC) ! total ice volume tot_volume = 0._dp do j = 1, nsn do i = 1, ewn if (model%geometry%thck(i,j) > model%numerics%thklim) then tot_volume = tot_volume + model%geometry%thck(i,j) endif enddo enddo tot_volume = tot_volume * model%numerics%dew * model%numerics%dns & * thk0 * len0**2 write(message,'(a25,e24.16)') 'Total ice volume (km^3) ', & tot_volume*1.0e-9_dp ! convert to km^3 call write_log(trim(message), type = GM_DIAGNOSTIC) ! total ice energy relative to T = 0 deg C tot_energy = 0._dp do j = 1, nsn do i = 1, ewn if (model%geometry%thck(i,j) > model%numerics%thklim) then thck = model%geometry%thck(i,j) do k = 1, upn-1 tot_energy = tot_energy & + thck * model%velowk%dups(k) * model%temper%temp(k,i,j) enddo endif enddo enddo tot_energy = tot_energy * model%numerics%dew * model%numerics%dns & * thk0 * len0**2 write(message,'(a25,e24.16)') 'Total ice energy (J) ', tot_energy call write_log(trim(message), type = GM_DIAGNOSTIC) ! total sum of volume * age tot_age = 0._dp do j = 1, nsn do i = 1, ewn if (model%geometry%thck(i,j) > model%numerics%thklim) then thck = model%geometry%thck(i,j) do k = 1,upn-1 tot_age = tot_age & + thck * model%velowk%dups(k) * model%geometry%age(k,i,j) enddo endif enddo enddo tot_age = tot_age * model%numerics%dew * model%numerics%dns & * thk0 * len0**2 * tim0/scyr ! mean thickness if (tot_area > eps) then mean_thck = tot_volume/tot_area else mean_thck = 0._dp endif write(message,'(a25,f24.16)') 'Mean thickness (m) ', mean_thck call write_log(trim(message), type = GM_DIAGNOSTIC) ! max thickness imax = 0 jmax = 0 max_thck = 0._dp do j = 1, nsn do i = 1, ewn if (model%geometry%thck(i,j) > max_thck) then max_thck = model%geometry%thck(i,j) imax = i jmax = j endif enddo enddo write(message,'(a25,f24.16,2i4)') 'Max thickness (m), i, j ', & max_thck*thk0, imax, jmax call write_log(trim(message), type = GM_DIAGNOSTIC) ! mean temperature if (tot_volume > eps) then mean_temp = tot_energy/tot_volume else mean_temp = 0._dp endif write(message,'(a25,f24.16)') 'Mean temperature (C) ', mean_temp call write_log(trim(message), type = GM_DIAGNOSTIC) ! min temperature min_temp = 999._dp do j = 1, nsn do i = 1, ewn do k = 1, upn-1 if (model%temper%temp(k,i,j) < min_temp) then min_temp = model%temper%temp(k,i,j) imin = i jmin = j kmin = k endif enddo enddo enddo write(message,'(a25,f24.16,3i4)') 'Min temperature, i, j, k ', & min_temp, imin, jmin, kmin call write_log(trim(message), type = GM_DIAGNOSTIC) ! mean age if (tot_age > eps) then mean_age = tot_age/tot_volume else mean_age = 0._dp endif write(message,'(a25,f24.16)') 'Mean ice age (yr) ', mean_age call write_log(trim(message), type = GM_DIAGNOSTIC) ! max surface speed ! Note: uvel and vvel not defined at i = ewn, j = nsn imax = 0 jmax = 0 max_spd_sfc = 0._dp k = 1 do j = 1, nsn-1 do i = 1, ewn-1 spd = sqrt(model%velocity%uvel(k,i,j)**2 & + model%velocity%vvel(k,i,j)**2) if (model%geometry%thck(i,j) > eps .and. spd > max_spd_sfc) then max_spd_sfc = spd imax = i jmax = j endif enddo enddo write(message,'(a25,f24.16,2i4)') 'Max sfc spd (m/yr), i, j ', & max_spd_sfc*vel0*scyr, imax, jmax call write_log(trim(message), type = GM_DIAGNOSTIC) ! max basal speed imax = 0 jmax = 0 max_spd_bas = 0._dp do j = 1, nsn-1 do i = 1, ewn-1 spd = sqrt(model%velocity%ubas(i,j)**2 & + model%velocity%vbas(i,j)**2) if (model%geometry%thck(i,j) > eps .and. spd > max_spd_bas) then max_spd_bas = spd imax = i jmax = j endif enddo enddo write(message,'(a25,f24.16,2i4)') 'Max base spd (m/yr), i, j', & max_spd_bas*vel0*scyr, imax, jmax call write_log(trim(message), type = GM_DIAGNOSTIC) ! local diagnostics if (present(idiag) .and. present(jdiag)) then if (idiag >= 1 .and. idiag <=model%general%ewn & .and. & jdiag >= 1 .and. jdiag <=model%general%nsn) then call write_log(' ') write(message,'(a30,2i4)') & 'Grid point diagnostics: i, j =', idiag, jdiag call write_log(trim(message), type = GM_DIAGNOSTIC) call write_log(' ') i = idiag j = jdiag write(message,'(a25,f24.16)') & 'Thickness (m) ',model%geometry%thck(i,j)*thk0 call write_log(trim(message), type = GM_DIAGNOSTIC) write(message,'(a25,f24.16)') & 'Sfc air temperature (C) ',model%climate%artm(i,j) call write_log(trim(message), type = GM_DIAGNOSTIC) write(message,'(a25,f24.16)') & 'Basal temperature (C) ',model%temper%temp(upn,i,j) call write_log(trim(message), type = GM_DIAGNOSTIC) spd = sqrt(model%velocity%ubas(i,j)**2 & + model%velocity%vbas(i,j)**2) * vel0*scyr write(message,'(a25,f24.16)') & 'Basal speed (m/yr) ', spd call write_log(trim(message), type = GM_DIAGNOSTIC) call write_log(' ') write(message,'(a52)') & 'Layer Speed (m/yr) Temperature (C) Age (yr) ' call write_log(trim(message), type = GM_DIAGNOSTIC) do k = 1, upn-1 spd = sqrt(model%velocity%uvel(k,i,j)**2 & + model%velocity%vvel(k,i,j)**2) * vel0*scyr write (message,'(i4,2f24.16)') & k, spd, model%temper%temp(k,i,j) ! age = model%geometry%age(k,i,j)*tim0/scyr ! uncomment when age is computed ! write (message,'(i4,3f24.16)') & ! k, spd, model%temper%temp(k,i,j), age call write_log(trim(message), type = GM_DIAGNOSTIC) enddo if (GLC_DEBUG) then do k = 1, upn-1 spd = sqrt(model%velocity%uvel(k,i,j)**2 & + model%velocity%vvel(k,i,j)**2) * vel0*scyr write (message,'(i4,Z24.20,Z24.20)') & k, spd, model%temper%temp(k,i,j) call write_log(trim(message), type = GM_DIAGNOSTIC) enddo endif endif ! idiag and jdiag in bounds endif ! idiag and jdiag present end subroutine glide_write_diag !============================================================== end module glide_diagnostics