diff --git a/models/cice-scm2/dart_cice_mod.f90 b/models/cice-scm2/dart_cice_mod.f90 deleted file mode 100644 index 5abe47e686..0000000000 --- a/models/cice-scm2/dart_cice_mod.f90 +++ /dev/null @@ -1,216 +0,0 @@ -! DART software - Copyright UCAR. This open source software is provided -! by UCAR, "as is", without charge, subject to all terms of use at -! http://www.image.ucar.edu/DAReS/DART/DART_download -! -! $Id$ - -module dart_cice_mod - -use types_mod, only : r8, rad2deg, PI, SECPERDAY, digits12 -use time_manager_mod, only : time_type, get_date, set_date, get_time, set_time, & - set_calendar_type, get_calendar_string, & - print_date, print_time, operator(==), operator(-) -use utilities_mod, only : get_unit, open_file, close_file, file_exist, & - register_module, error_handler, & - find_namelist_in_file, check_namelist_read, & - E_ERR, E_MSG, find_textfile_dims - -use netcdf_utilities_mod, only : nc_check - - -use typesizes -use netcdf - -implicit none -private - -public :: set_model_time_step,get_horiz_grid_dims, & - get_ncat_dim, read_horiz_grid - -character(len=*), parameter :: source = "$URL$" -character(len=*), parameter :: revision = "$Revision$" -character(len=*), parameter :: revdate = "$Date$" - -character(len=512) :: msgstring -logical, save :: module_initialized = .false. - -character(len=256) :: ic_filename = 'cice.r.nc' - -contains - -subroutine initialize_module - -integer :: iunit, io - -! Read calendar information -! In 'restart' mode, this is primarily the calendar type and 'stop' -! information. The time attributes of the restart file override -! the namelist time information. - -! FIXME : Real observations are always GREGORIAN dates ... -! but stomping on that here gets in the way of running -! a perfect_model experiment for pre-1601 AD cases. -call set_calendar_type('gregorian') - -! Make sure we have a cice restart file (for grid dims) -if ( .not. file_exist(ic_filename) ) then - msgstring = 'dart_cice_mod: '//trim(ic_filename)//' not found' - call error_handler(E_ERR,'initialize_module', & - msgstring, source, revision, revdate) -endif - -module_initialized = .true. - -! Print module information to log file and stdout. -call register_module(source, revision, revdate) - -end subroutine initialize_module -!!!!!!!!!!!!!!!! -function set_model_time_step() - -! the initialize_module ensures that the cice namelists are read. -! The restart times in the cice_in&restart_nml are used to define -! appropriate assimilation timesteps. -! -type(time_type) :: set_model_time_step - -if ( .not. module_initialized ) call initialize_module - -! Check the 'restart_option' and 'restart_n' to determine -! when we can stop the model -! CMB not sure if nday is actually different than ndays, no matter here though -!if ( (trim(restart_option) == 'ndays') .or. (trim(restart_option) == 'nday' ) ) then -! set_model_time_step = set_time(0, restart_n) ! (seconds, days) -!else if ( trim(restart_option) == 'nyears' ) then - ! FIXME ... CMB I guess we ignore it and make the freq 1 day anyway? - set_model_time_step = set_time(0, 1) ! (seconds, days) -!else -! call error_handler(E_ERR,'set_model_time_step', & -! 'restart_option must be ndays or nday', source, revision, revdate) -!endif - -end function set_model_time_step -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine get_horiz_grid_dims(Nx) - -! -! Read the lon, lat grid size from the restart netcdf file. -! The actual grid file is a binary file with no header information. -! -! The file name comes from module storage ... namelist. - -integer, intent(out) :: Nx ! Number of Longitudes - -integer :: grid_id, dimid, nc_rc - -if ( .not. module_initialized ) call initialize_module - -call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), & - 'get_horiz_grid_dims','open '//trim(ic_filename)) - -! Longitudes : get dimid for 'ni' or 'nlon', and then get value -nc_rc = nf90_inq_dimid(grid_id, 'ni', dimid) -if (nc_rc /= nf90_noerr) then - msgstring = "unable to find either 'ni' or 'nlon' in file "//trim(ic_filename) - call error_handler(E_ERR, 'get_horiz_grid_dims', msgstring, & - source,revision,revdate) -endif - -call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Nx), & - 'get_horiz_grid_dims','inquire_dimension ni '//trim(ic_filename)) - -call nc_check(nf90_close(grid_id), & - 'get_horiz_grid_dims','close '//trim(ic_filename) ) - -end subroutine get_horiz_grid_dims -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine get_ncat_dim(Ncat) - -! -! Read the ncat size from the restart netcdf file. - -integer, intent(out) :: Ncat ! Number of categories in ice-thick dist - -integer :: grid_id, dimid, nc_rc - -if ( .not. module_initialized ) call initialize_module - -call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), & - 'get_ncat_dim','open '//trim(ic_filename)) - -! ncat : get dimid for 'ncat' and then get value -nc_rc = nf90_inq_dimid(grid_id, 'ncat', dimid) -if (nc_rc /= nf90_noerr) then - nc_rc = nf90_inq_dimid(grid_id, 'Ncat', dimid) - if (nc_rc /= nf90_noerr) then - msgstring = "unable to find either 'ncat' or 'Ncat' in file "//trim(ic_filename) - call error_handler(E_ERR, 'get_horiz_grid_dims', msgstring, & - source,revision,revdate) - endif -endif - -call nc_check(nf90_inquire_dimension(grid_id, dimid, len=Ncat), & - 'get_ncat_dim','inquire_dimension ni '//trim(ic_filename)) - -! tidy up - -call nc_check(nf90_close(grid_id), & - 'get_ncat_dim','close '//trim(ic_filename) ) - -end subroutine get_ncat_dim -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine read_horiz_grid(nx, TLAT, TLON) - -integer, intent(in) :: nx -real(r8), dimension(nx), intent(out) :: TLAT, TLON - -integer :: grid_id, reclength,VarId,status - -if ( .not. module_initialized ) call initialize_module - -! Check to see that the file exists. - -if ( .not. file_exist(ic_filename) ) then - msgstring = 'cice grid '//trim(ic_filename)//' not found' - call error_handler(E_ERR,'read_horiz_grid', & - msgstring, source, revision, revdate) -endif - -! Open it and read them in the EXPECTED order. -! Actually, we only need the first two, so I'm skipping the rest. - -call nc_check(nf90_open(trim(ic_filename), nf90_nowrite, grid_id), & - 'read_horiz_grid','open '//trim(ic_filename)) -! Latitude -call nc_check(nf90_inq_varid(grid_id, 'tlat', VarId), & - 'read_horiz_grid','inquiring tlat from '//trim(ic_filename)) -call nc_check(nf90_get_var(grid_id, VarId, TLAT, & - start=(/1/), & - count=(/nx/)), & -'read_horiz_grid','getting tlat from '//trim(ic_filename)) -!Longitude -call nc_check(nf90_inq_varid(grid_id, 'tlon', VarId), & -'read_horiz_grid','inquiring tlon from '//trim(ic_filename)) -call nc_check(nf90_get_var(grid_id, VarId, TLON, & - start=(/1/), & - count=(/nx/)), & - 'read_horiz_grid','getting tlon from '//trim(ic_filename)) - -call nc_check(nf90_close(grid_id), & - 'read_horiz_grid','close '//trim(ic_filename) ) - -TLAT = TLAT * rad2deg -TLON = TLON * rad2deg - -! ensure [0,360) [-90,90] - -where (TLON < 0.0_r8) TLON = TLON + 360.0_r8 -where (TLON > 360.0_r8) TLON = TLON - 360.0_r8 - -where (TLAT < -90.0_r8) TLAT = -90.0_r8 -where (TLAT > 90.0_r8) TLAT = 90.0_r8 - -end subroutine read_horiz_grid - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module dart_cice_mod diff --git a/models/cice-scm2/dart_to_cice.f90 b/models/cice-scm2/dart_to_cice.f90 deleted file mode 100644 index 3882e2fd7c..0000000000 --- a/models/cice-scm2/dart_to_cice.f90 +++ /dev/null @@ -1,578 +0,0 @@ -! DART software - Copyright UCAR. This open source software is provided -! by UCAR, "as is", without charge, subject to all terms of use at -! http://www.image.ucar.edu/DAReS/DART/DART_download -! -! $Id$ - -program dart_to_cice - -!---------------------------------------------------------------------- -! purpose: implement a 'partition function' to modify the cice state -! to be consistent with the states from assimilation -! -! method: Read in restart (restart with prior) and out restart (restart -! with posterior) written by DART after filter. -! -! author: C Bitz June 2016 -!---------------------------------------------------------------------- - -use types_mod, only : r8 -use utilities_mod, only : initialize_utilities, finalize_utilities, & - find_namelist_in_file, check_namelist_read, & - file_exist, error_handler, E_ERR, E_MSG, to_upper -use netcdf_utilities_mod, only : nc_check -use netcdf - - -implicit none - -! version controlled file description for error handling, do not edit -character(len=*), parameter :: source = & - "$URL$" -character(len=*), parameter :: revision = "$Revision$" -character(len=*), parameter :: revdate = "$Date$" - -!------------------------------------------------------------------ - -character(len=256) :: dart_to_cice_input_file = 'dart_restart.nc' -character(len=256) :: original_cice_input_file = 'cice_restart.nc' -character(len=256) :: previous_cice_input_file = 'pre_restart.nc' -character(len=128) :: balance_method = 'simple_squeeze' -character(len=15) :: r_snw_name = 'r_snw' -integer :: gridpt_oi = 3 - -namelist /dart_to_cice_nml/ dart_to_cice_input_file, & - original_cice_input_file, & - previous_cice_input_file, & - balance_method, & - r_snw_name, & - gridpt_oi - -character(len=512) :: string1, string2, msgstring -character(len=15) :: varname -character(len=128) :: method - -integer :: Nx -integer :: Ncat ! number of categories in ice-thickness dist -integer, parameter :: Nilyr = 8 ! number of layers in ice, hardwired -integer, parameter :: Nslyr = 3 ! number of layers in snow, hardwired - -real(r8), allocatable :: aicen_original(:) -real(r8), allocatable :: vicen_original(:) -real(r8), allocatable :: vsnon_original(:) -!real(r8), allocatable :: aice_original(:,:) -!real(r8), allocatable :: hicen_original(:) -!real(r8), allocatable :: hsnon_original(:) -logical :: sst_present = .true. -logical :: sst_org_present = .true. - -real(r8) :: sst,sst_original -real(r8), allocatable :: aicen(:) -real(r8), allocatable :: vicen(:) -real(r8), allocatable :: vsnon(:) -real(r8), allocatable :: Tsfcn(:) -real(r8), allocatable :: qice(:,:) -real(r8), allocatable :: sice(:,:) -real(r8), allocatable :: qsno(:,:) - -character (len=3) :: nchar -integer :: iunit,io,ncid,dimid,l,n,VarID -real(r8) :: aice,aice_temp -real(r8) :: vice,vice_temp -real(r8) :: vsno,vsno_temp -real(r8), parameter :: Tsmelt = 0._r8 -real(r8), parameter :: c1 = 1.0_r8 -real(r8), parameter :: & - phi_init = 0.75_r8, & - dSin0_frazil = 3.0_r8 -real(r8), parameter :: sss = 34.7_r8 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -real(r8) :: squeeze,cc1,cc2,cc3,x1,Si0new,Ti,qsno_hold,qi0new -real(r8), allocatable :: hin_max(:) -real(r8), allocatable :: hcat_midpoint(:) - -call initialize_utilities(progname='dart_to_cice') - -call find_namelist_in_file("input.nml", "dart_to_cice_nml", iunit) -read(iunit, nml = dart_to_cice_nml, iostat = io) -call check_namelist_read(iunit, io, "dart_to_cice_nml") - -method = balance_method -call to_upper(method) - -! check on namelist stuff, and whether files exist -write(string1,*) 'converting DART output file "'// & - &trim(dart_to_cice_input_file)//'" to one CICE will like' -write(string2,*) 'using the "'//trim(balance_method)//'" method.' -call error_handler(E_MSG,'dart_to_cice',string1,text2=string2) - -if ( .not. file_exist(dart_to_cice_input_file) ) then - write(string1,*) 'cannot open "', trim(dart_to_cice_input_file),'" for updating.' - call error_handler(E_ERR,'dart_to_cice:filename not found ',trim(dart_to_cice_input_file)) -endif - -if ( .not. file_exist(original_cice_input_file) ) then - write(string1,*) 'cannot open "', trim(original_cice_input_file),'" for reading.' - call error_handler(E_ERR,'dart_to_cice:filename not found ',trim(original_cice_input_file)) -endif - - -call nc_check( nf90_open(trim(original_cice_input_file), NF90_NOWRITE, ncid), & - 'dart_to_cice', 'open "'//trim(original_cice_input_file)//'"') - -call nc_check(nf90_inq_dimid(ncid,"ncat",dimid), & - 'dart_to_cice', 'inquire ncat dimid from "'//trim(original_cice_input_file)//'"') -call nc_check(nf90_inquire_dimension(ncid,dimid,len=Ncat), & - 'dart_to_cice', 'inquire ncat from "'//trim(original_cice_input_file)//'"') -call nc_check(nf90_inq_dimid(ncid,"ni",dimid), & - 'dart_to_cice', 'inquire ni dimid from "'//trim(original_cice_input_file)//'"') -call nc_check(nf90_inquire_dimension(ncid,dimid,len=Nx),& - 'dart_to_cice', 'inquire ni from "'//trim(original_cice_input_file)//'"') - -allocate(aicen_original(NCAT),vicen_original(NCAT),vsnon_original(NCAT),Tsfcn(NCAT),qice(Nilyr,NCAT),sice(Nilyr,NCAT),qsno(Nslyr,NCAT)) -call get_variable(ncid,'aicen',aicen_original,original_cice_input_file,gridpt_oi,Ncat) -call get_variable(ncid,'vicen',vicen_original,original_cice_input_file,gridpt_oi,Ncat) -call get_variable(ncid,'vsnon',vsnon_original,original_cice_input_file,gridpt_oi,Ncat) -call get_variable(ncid,'Tsfcn',Tsfcn,dart_to_cice_input_file,gridpt_oi,Ncat) -call get_variable1d(ncid,'sst',sst_original,dart_to_cice_input_file,gridpt_oi,sst_org_present) -do l=1, Nilyr - write(nchar,'(i3.3)') l - call get_variable(ncid,'qice'//trim(nchar),qice(l,:),dart_to_cice_input_file,gridpt_oi,Ncat) - call get_variable(ncid,'sice'//trim(nchar),sice(l,:),dart_to_cice_input_file,gridpt_oi,Ncat) -enddo -do l=1, Nslyr - write(nchar,'(i3.3)') l - call get_variable(ncid,'qsno'//trim(nchar),qsno(l,:),dart_to_cice_input_file,gridpt_oi,Ncat) -enddo -call nc_check(nf90_close(ncid),'dart_to_cice', 'close '//trim(original_cice_input_file)) -!!!!!!!!! -call nc_check( nf90_open(trim(dart_to_cice_input_file), NF90_NOWRITE, ncid), & - 'dart_to_cice', 'open "'//trim(dart_to_cice_input_file)//'"') -allocate(aicen(NCAT),vicen(NCAT),vsnon(NCAT)) -call get_variable(ncid,'aicen',aicen,dart_to_cice_input_file,gridpt_oi,Ncat) -call get_variable(ncid,'vicen',vicen,dart_to_cice_input_file,gridpt_oi,Ncat) -call get_variable(ncid,'vsnon',vsnon,dart_to_cice_input_file,gridpt_oi,Ncat) -call get_variable1d(ncid,'sst',sst,dart_to_cice_input_file,gridpt_oi,sst_present) -call nc_check(nf90_close(ncid),'dart_to_cice', 'close '//trim(dart_to_cice_input_file)) -!!!!!!!!!!!!!!!!!!!!!!!!! -qice = min(0.0_r8,qice) -sice = max(0.0_r8,sice) -qsno = min(0.0_r8,qsno) -aicen = min(1.0_r8,aicen) -Tsfcn = min(Tsmelt,Tsfcn) -!!!!!! -aice = sum(aicen) -vice = sum(vicen) -vsno = sum(vsnon) -!!!!!! -aicen = max(0.0_r8,aicen) -vicen = max(0.0_r8,vicen) -vsnon = max(0.0_r8,vsnon) -!!!!! -aice_temp = sum(aicen) -vice_temp = sum(vicen) -vsno_temp = sum(vsnon) -!!!!! -if (aice<0.0_r8) then - aicen(:) = 0.0_r8 - vicen(:) = 0.0_r8 - vsnon(:) = 0.0_r8 -endif -!!!!! -do n=1,NCAT - if (aice_temp > 0._r8 .and. aice>0._r8) then - aicen(n) = aicen(n) - (aice_temp-aice)*aicen(n)/aice_temp - endif - if (vice_temp > 0._r8 .and. vice>0._r8) then - vicen(n) = vicen(n) - (vice_temp-vice)*vicen(n)/vice_temp - endif - if (vsno_temp > 0._r8 .and. vsno > 0._r8) then - vsnon(n) = vsnon(n) - (vsno_temp-vsno)*vsnon(n)/vsno_temp - endif -enddo -!!!! -if (aice>1.0_r8) then - squeeze = 1.0_r8/aice - aicen(:) = aicen(:)*squeeze -endif -!!!!!! -if (sst_present) then - if (aice == 0.0_r8) sst = 0.0_r8 -endif -where(aicen==-999) aicen = 0.0_r8 -!!!!!! -cc1 = 3._r8/real(Ncat,kind=r8) -cc2 = 15.0_r8*cc1 -cc3 = 3._r8 -allocate( hin_max(0:Ncat) ) -allocate( hcat_midpoint(Ncat) ) -hin_max(0) = 0._r8 -do n = 1, NCAT - x1 = real(n-1,kind=r8) / real(Ncat,kind=r8) - hin_max(n) = hin_max(n-1) & - + cc1 + cc2*(c1 + tanh(cc3*(x1-c1))) - hcat_midpoint(n)=0.5_r8*(hin_max(n-1)+hin_max(n)) -enddo -!!!!!!! -do n=1,NCAT - if (aicen(n) > 0.0_r8 .and. aicen_original(n) > 0.0_r8) then - if (vicen(n) == 0.0_r8) then - vicen(n) = aicen(n)*hcat_midpoint(n) - endif - endif - if (aicen(n) == 0.0_r8 .and. aicen_original(n) > 0.0_r8) then - vicen(n) = 0.0_r8 - qice(:,n) = 0.0_r8 - sice(:,n) = 0.0_r8 - qsno(:,n) = 0.0_r8 - vsnon(n) = 0.0_r8 - Tsfcn(n) = -1.8_r8 - else if (aicen(n)>0.0_r8 .and. aicen_original(n) == 0.0_r8) then - if (vicen(n) == 0.0_r8) vicen(n) = aicen(n) * hcat_midpoint(n) - Si0new = sss - dSin0_frazil - sice(:,n) = Si0new - Ti = min(liquidus_temperature_mush(Si0new/phi_init), -0.1_r8) - qi0new = enthalpy_mush(Ti, Si0new) - qice(:,n) = qi0new - if (vsnon(n) == 0.0_r8 .and. vsnon_original(n) > 0.0_r8) then - qsno(:,n) = 0.0_r8 - else if (vsnon(n) > 0.0_r8 .and. vsnon_original(n) == 0.0_r8) then - qsno_hold = snow_enthaply(Ti) - qsno(:,n) = qsno_hold - endif - Tsfcn(n) = Ti - endif - if (aicen(n) == 0.0_r8) then - vicen(n) = 0.0_r8 - vsnon(n) = 0.0_r8 - endif -enddo -!!!!!!!! -call nc_check( nf90_open(trim(original_cice_input_file), NF90_WRITE, ncid), & - 'dart_to_cice', 'open "'//trim(original_cice_input_file)//'"') -varname='aicen' -io = nf90_inq_varid(ncid, trim(varname), VarID) -call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) -io = nf90_put_var(ncid, VarID, aicen,start=(/gridpt_oi,1/),count=(/1,NCAT/)) -call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -!!!! -varname='vicen' -io = nf90_inq_varid(ncid, trim(varname), VarID) -call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) -io = nf90_put_var(ncid, VarID, vicen,start=(/gridpt_oi,1/),count=(/1,NCAT/)) -call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -!!!! -varname='vsnon' -io = nf90_inq_varid(ncid, trim(varname), VarID) -call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) -io = nf90_put_var(ncid, VarID, vsnon,start=(/gridpt_oi,1/),count=(/1,NCAT/)) -call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -!!!! -varname='Tsfcn' -io = nf90_inq_varid(ncid, trim(varname), VarID) -call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) -io = nf90_put_var(ncid, VarID, Tsfcn,start=(/gridpt_oi,1/),count=(/1,NCAT/)) -call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -!!!!! -if (sst_present) then - varname='sst' - io = nf90_inq_varid(ncid, trim(varname), VarID) - call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) - io = nf90_put_var(ncid, VarID, sst,start=(/gridpt_oi/))!,count=(/1/)) - call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -endif -!!!!! -do l=1, Nilyr - write(nchar,'(i3.3)') l - varname='qice'//trim(nchar) - io = nf90_inq_varid(ncid, trim(varname), VarID) - call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) - io = nf90_put_var(ncid, VarID, qice(l,:),start=(/gridpt_oi,1/),count=(/1,NCAT/)) - call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) - !!!!!!!!!! - varname='sice'//trim(nchar) - io = nf90_inq_varid(ncid, trim(varname), VarID) - call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) - io = nf90_put_var(ncid, VarID, sice(l,:),start=(/gridpt_oi,1/),count=(/1,NCAT/)) - call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -enddo -!!!! -do l=1, Nslyr - write(nchar,'(i3.3)') l - varname='qsno'//trim(nchar) - io = nf90_inq_varid(ncid, trim(varname), VarID) - call nc_check(io, 'dart_to_cice', & - 'inq_varid '//trim(varname)//' '//trim(original_cice_input_file)) - io = nf90_put_var(ncid, VarID, qsno(l,:),start=(/gridpt_oi,1/),count=(/1,NCAT/)) - call nc_check(io, 'dart_to_cice', & - 'put_var '//trim(varname)//' '//trim(original_cice_input_file)) -enddo - -call nc_check(nf90_close(ncid),'dart_to_cice', 'close '//trim(original_cice_input_file)) - -deallocate( aicen, vicen, vsnon, Tsfcn) -deallocate( qice, sice, qsno ) - - -call finalize_utilities('dart_to_cice') - - -contains - -subroutine get_variable(ncid,varname,var,filename,space_index,ncat) -integer, intent(in) :: ncid,ncat -character(len=*), intent(in) :: varname -real(r8), intent(out) :: var(ncat) -character(len=*), intent(in) :: filename -integer, intent(in) :: space_index - -integer :: VarID, ndims, dimIDs -real(r8) :: holder(4,ncat) - -write(6,*) 'Getting data for ',trim(varname) - -io = nf90_inq_varid(ncid, trim(varname), VarID) -call nc_check(io, 'dart_to_cice', 'inq_varid '//trim(msgstring)) - -call nc_check(nf90_get_var(ncid, VarID, holder), 'dart_to_cice', & - 'get_var '//trim(msgstring)) - - -var(:) = holder(gridpt_oi,:) - -end subroutine get_variable -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine get_variable1d(ncid,varname,var,filename,space_index,var_present) -integer, intent(in) :: ncid -character(len=*), intent(in) :: varname -real(r8), intent(out) :: var -character(len=*), intent(in) :: filename -integer, intent(in) :: space_index -logical, intent(inout) :: var_present - -integer :: VarID, ndims, dimIDs -real(r8) :: holder(4) - -write(6,*) 'Getting data for ',trim(varname) - -io = nf90_inq_varid(ncid, trim(varname), VarID) -if(io /= nf90_NoErr) then - write(6,*) "No netcdf ID for ",trim(varname) - var_present = .false. - return -endif -call nc_check(io, 'dart_to_cice', 'inq_varid '//trim(msgstring)) - -call nc_check(nf90_get_var(ncid, VarID, holder), 'dart_to_cice', & - 'get_var '//trim(msgstring)) - - -var = holder(gridpt_oi) - -end subroutine get_variable1d -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function enthalpy_mush(zTin, zSin) result(zqin) - - ! enthalpy of mush from mush temperature and bulk salinity - - real(r8), intent(in) :: & - zTin, & ! ice layer temperature (C) - zSin ! ice layer bulk salinity (ppt) - - real(r8) :: & - zqin ! ice layer enthalpy (J m-3) - - real(r8) :: & - phi ! ice liquid fraction - -! from shr_const_mod.F90 - real(r8),parameter :: SHR_CONST_CPSW = 3.996e3_R8 ! specific heat of sea water ~ J/kg/K - real(R8),parameter :: SHR_CONST_CPICE = 2.11727e3_R8 ! specific heat of fresh ice ~ J/kg/K - real(R8),parameter :: SHR_CONST_RHOSW = 1.026e3_R8 ! density of sea water ~ kg/m^3 - real(R8),parameter :: SHR_CONST_RHOICE= 0.917e3_R8 ! density of ice ~ kg/m^3 - real(R8),parameter :: SHR_CONST_LATICE= 3.337e5_R8 ! latent heat of fusion ~ J/kg - - -! from cice/src/drivers/cesm/ice_constants.F90 - real(r8) :: cp_ocn, cp_ice, rhoi, rhow, Lfresh - - cp_ice = SHR_CONST_CPICE ! specific heat of fresh ice (J/kg/K) - cp_ocn = SHR_CONST_CPSW ! specific heat of ocn (J/kg/K) - rhoi = SHR_CONST_RHOICE ! density of ice (kg/m^3) - rhow = SHR_CONST_RHOSW ! density of seawater (kg/m^3) - Lfresh = SHR_CONST_LATICE ! latent heat of melting of fresh ice (J/kg) - - phi = liquid_fraction(zTin, zSin) - - zqin = phi * (cp_ocn * rhow - cp_ice * rhoi) * zTin + & - rhoi * cp_ice * zTin - (1._r8 - phi) * rhoi * Lfresh - - end function enthalpy_mush -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function liquid_fraction(zTin, zSin) result(phi) - - ! liquid fraction of mush from mush temperature and bulk salinity - - real(r8), intent(in) :: & - zTin, & ! ice layer temperature (C) - zSin ! ice layer bulk salinity (ppt) - - real(r8) :: & - phi , & ! liquid fraction - Sbr ! brine salinity (ppt) - - real (r8), parameter :: puny = 1.0e-11_r8 ! cice/src/drivers/cesm/ice_constants.F90 - - Sbr = max(liquidus_brine_salinity_mush(zTin),puny) - phi = zSin / max(Sbr, zSin) - - end function liquid_fraction -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function snow_enthaply(Ti) result(qsno) - real(r8), intent(in) :: Ti - - real(r8),parameter :: rhos = 330.0_r8, & - Lfresh = 2.835e6_r8 - 2.501e6_r8, & - cp_ice = 2106._r8 - real(r8) :: qsno - - qsno = -rhos*(Lfresh - cp_ice*min(0.0_r8,Ti)) - end function snow_enthaply -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function liquidus_brine_salinity_mush(zTin) result(Sbr) - - ! liquidus relation: equilibrium brine salinity as function of temperature - ! based on empirical data from Assur (1958) - - real(r8), intent(in) :: & - zTin ! ice layer temperature (C) - - real(r8) :: & - Sbr ! ice brine salinity (ppt) - - real(r8) :: & - t_high , & ! mask for high temperature liquidus region - lsubzero ! mask for sub-zero temperatures - - !constant numbers from ice_constants.F90 - real(r8), parameter :: & - c1 = 1.0_r8 , & - c1000 = 1000_r8 - - ! liquidus relation - higher temperature region - real(r8), parameter :: & - az1_liq = -18.48_r8 ,& - bz1_liq = 0.0_r8 - - ! liquidus relation - lower temperature region - real(r8), parameter :: & - az2_liq = -10.3085_r8, & - bz2_liq = 62.4_r8 - - ! liquidus break - real(r8), parameter :: & - Tb_liq = -7.6362968855167352_r8 - - ! basic liquidus relation constants - real(r8), parameter :: & - az1p_liq = az1_liq / c1000, & - bz1p_liq = bz1_liq / c1000, & - az2p_liq = az2_liq / c1000, & - bz2p_liq = bz2_liq / c1000 - - ! temperature to brine salinity - real(r8), parameter :: & - J1_liq = bz1_liq / az1_liq , & - K1_liq = c1 / c1000 , & - L1_liq = (c1 + bz1p_liq) / az1_liq , & - J2_liq = bz2_liq / az2_liq , & - K2_liq = c1 / c1000 , & - L2_liq = (c1 + bz2p_liq) / az2_liq - - t_high = merge(1._r8, 0._r8, (zTin > Tb_liq)) - lsubzero = merge(1._r8, 0._r8, (zTin <= 1._r8)) - - Sbr = ((zTin + J1_liq) / (K1_liq * zTin + L1_liq)) * t_high + & - ((zTin + J2_liq) / (K2_liq * zTin + L2_liq)) * (1._r8 - t_high) - - Sbr = Sbr * lsubzero - - end function liquidus_brine_salinity_mush -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function liquidus_temperature_mush(Sbr) result(zTin) - - ! liquidus relation: equilibrium temperature as function of brine salinity - ! based on empirical data from Assur (1958) - - real(r8), intent(in) :: & - Sbr ! ice brine salinity (ppt) - - real(r8) :: & - zTin ! ice layer temperature (C) - - real(r8) :: & - t_high ! mask for high temperature liquidus region - - ! liquidus break - real(r8), parameter :: & - Sb_liq = 123.66702800276086_r8 ! salinity of liquidus break - - ! constant numbers from ice_constants.F90 - real(r8), parameter :: & - c1 = 1.0_r8 , & - c1000 = 1000_r8 - - ! liquidus relation - higher temperature region - real(r8), parameter :: & - az1_liq = -18.48_r8 ,& - bz1_liq = 0.0_r8 - - ! liquidus relation - lower temperature region - real(r8), parameter :: & - az2_liq = -10.3085_r8, & - bz2_liq = 62.4_r8 - - ! basic liquidus relation constants - real(r8), parameter :: & - az1p_liq = az1_liq / c1000, & - bz1p_liq = bz1_liq / c1000, & - az2p_liq = az2_liq / c1000, & - bz2p_liq = bz2_liq / c1000 - - ! brine salinity to temperature - real(r8), parameter :: & - M1_liq = az1_liq , & - N1_liq = -az1p_liq , & - O1_liq = -bz1_liq / az1_liq , & - M2_liq = az2_liq , & - N2_liq = -az2p_liq , & - O2_liq = -bz2_liq / az2_liq - - t_high = merge(1._r8, 0._r8, (Sbr <= Sb_liq)) - - zTin = ((Sbr / (M1_liq + N1_liq * Sbr)) + O1_liq) * t_high + & - ((Sbr / (M2_liq + N2_liq * Sbr)) + O2_liq) * (1._r8 - t_high) - - end function liquidus_temperature_mush -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end program dart_to_cice - -! -! $URL$ -! $Id$ -! $Revision$ -! $Date$ diff --git a/models/cice-scm2/model_mod.f90 b/models/cice-scm2/model_mod.f90 deleted file mode 100644 index a2c1b617ed..0000000000 --- a/models/cice-scm2/model_mod.f90 +++ /dev/null @@ -1,1071 +0,0 @@ -! DART software - Copyright UCAR. This open source software is provided -! by UCAR, "as is", without charge, subject to all terms of use at -! http://www.image.ucar.edu/DAReS/DART/DART_download -! -! $Id$ - -module model_mod - -! This is a template showing the interfaces required for a model to be compliant -! with the DART data assimilation infrastructure. The public interfaces listed -! must all be supported with the argument lists as indicated. Many of the interfaces -! are not required for minimal implementation (see the discussion of each -! interface and look for NULL INTERFACE). - -! Modules that are absolutely required for use are listed -use types_mod, only : i4, r8, i8, MISSING_R8, metadatalength -use time_manager_mod, only : time_type, set_time, set_time_missing,set_calendar_type,get_time, & - set_date, get_date -use location_mod, only : location_type, get_close_type, & - get_close_obs, get_dist,& - convert_vertical_obs, convert_vertical_state, & - set_location, set_location_missing,VERTISLEVEL, & - get_location, & - loc_get_close_state => get_close_state -use utilities_mod, only : register_module, error_handler, & - E_ERR, E_MSG, logfileunit, & - nmlfileunit, do_output, do_nml_file, do_nml_term, & - find_namelist_in_file, check_namelist_read,to_upper, & - file_exist -use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & - nc_add_global_creation_time, & - nc_begin_define_mode, nc_end_define_mode, & - nc_check -use state_structure_mod, only : add_domain, get_domain_size -use ensemble_manager_mod, only : ensemble_type -use distributed_state_mod, only : get_state -use default_model_mod, only : pert_model_copies, nc_write_model_vars, init_conditions, & - init_time, adv_1step -use dart_cice_mod, only : set_model_time_step,get_horiz_grid_dims, & - get_ncat_dim, read_horiz_grid -use state_structure_mod, only : state_structure_info,get_index_start, get_num_variables, & - get_dart_vector_index, get_model_variable_indices -use obs_kind_mod, only : QTY_SEAICE_AGREG_CONCENTR , & - QTY_SEAICE_AGREG_VOLUME , & - QTY_SEAICE_AGREG_SNOWVOLUME, & - QTY_SEAICE_AGREG_THICKNESS , & - QTY_SEAICE_AGREG_SNOWDEPTH , & - QTY_SEAICE_CATEGORY , & - QTY_U_SEAICE_COMPONENT , & - QTY_V_SEAICE_COMPONENT , & - QTY_SEAICE_ALBEDODIRVIZ , & - QTY_SEAICE_ALBEDODIRNIR , & - QTY_SEAICE_ALBEDOINDVIZ , & - QTY_SEAICE_ALBEDOINDNIR , & - QTY_SEAICE_CONCENTR , & - QTY_SEAICE_VOLUME , & - QTY_SEAICE_SNOWVOLUME , & - QTY_SEAICE_SURFACETEMP , & - QTY_SEAICE_FIRSTYEARAREA , & - QTY_SEAICE_ICEAGE , & - QTY_SEAICE_LEVELAREA , & - QTY_SEAICE_LEVELVOLUME , & - QTY_SEAICE_MELTPONDAREA , & - QTY_SEAICE_MELTPONDDEPTH , & - QTY_SEAICE_MELTPONDLID , & - QTY_SEAICE_MELTPONDSNOW , & - QTY_SEAICE_SALINITY001 , & - QTY_SEAICE_SALINITY002 , & - QTY_SEAICE_SALINITY003 , & - QTY_SEAICE_SALINITY004 , & - QTY_SEAICE_SALINITY005 , & - QTY_SEAICE_SALINITY006 , & - QTY_SEAICE_SALINITY007 , & - QTY_SEAICE_SALINITY008 , & - QTY_SEAICE_ICEENTHALPY001 , & - QTY_SEAICE_ICEENTHALPY002 , & - QTY_SEAICE_ICEENTHALPY003 , & - QTY_SEAICE_ICEENTHALPY004 , & - QTY_SEAICE_ICEENTHALPY005 , & - QTY_SEAICE_ICEENTHALPY006 , & - QTY_SEAICE_ICEENTHALPY007 , & - QTY_SEAICE_ICEENTHALPY008 , & - QTY_SEAICE_SNOWENTHALPY001 , & - QTY_SEAICE_SNOWENTHALPY002 , & - QTY_SEAICE_SNOWENTHALPY003 , & - QTY_DRY_LAND , & - QTY_SOM_TEMPERATURE , & - QTY_SEAICE_FY , & - QTY_SEAICE_AGREG_FY , & - QTY_SEAICE_AGREG_SURFACETEMP,& - get_index_for_quantity , & - get_name_for_quantity - -use netcdf - -implicit none -private - -! required by DART code - will be called from filter and other -! DART executables. interfaces to these routines are fixed and -! cannot be changed in any way. -public :: get_model_size, & - adv_1step, & - get_state_meta_data, & - model_interpolate, & - shortest_time_between_assimilations, & - end_model, & - static_init_model, & - nc_write_model_atts, & - init_time, & - init_conditions, & - check_sfctemp_var - -! public but in another module -public :: nc_write_model_vars, & - pert_model_copies, & - get_close_obs, & - get_close_state, & - convert_vertical_obs, & - convert_vertical_state, & - read_model_time, & - write_model_time - - -! version controlled file description for error handling, do not edit -character(len=256), parameter :: source = & - "$URL$" -character(len=32 ), parameter :: revision = "$Revision$" -character(len=128), parameter :: revdate = "$Date$" -character(len=512) :: string1 -character(len=512) :: string2 -character(len=512) :: string3 - -type(location_type), allocatable :: state_loc(:) ! state locations, compute once and store for speed - -type(time_type) :: assimilation_time_step - -! DART state vector contents are specified in the input.nml:&model_nml namelist. -integer, parameter :: max_state_variables = 10 -integer, parameter :: num_state_table_columns = 3 -character(len=NF90_MAX_NAME) :: variable_table( max_state_variables, num_state_table_columns ) -integer :: state_kinds_list( max_state_variables ) -logical :: update_var_list( max_state_variables ) - -integer, parameter :: VAR_NAME_INDEX = 1 -integer, parameter :: VAR_QTY_INDEX = 2 -integer, parameter :: VAR_UPDATE_INDEX = 3 - -! EXAMPLE: perhaps a namelist here for anything you want to/can set at runtime. -! this is optional! only add things which can be changed at runtime. -integer :: model_size -integer :: assimilation_period_days = 0 -integer :: assimilation_period_seconds = 3600 - -real(r8) :: model_perturbation_amplitude = 0.01 - -character(len=metadatalength) :: model_state_variables(max_state_variables * num_state_table_columns ) = ' ' -integer :: debug = 100 -integer :: grid_oi = 3 -logical, save :: module_initialized = .false. - -real(r8), allocatable :: TLAT(:), TLON(:) - -type(time_type) :: model_time, model_timestep - -integer :: Nx=-1 -integer :: Ncat=-1 -integer :: domain_id,nfields -! uncomment this, the namelist related items in the 'use utilities' section above, -! and the namelist related items below in static_init_model() to enable the -! run-time namelist settings. -!namelist /model_nml/ model_size, assimilation_time_step_days, assimilation_time_step_seconds - -namelist /model_nml/ & - assimilation_period_days, & ! for now, this is the timestep - assimilation_period_seconds, & - model_perturbation_amplitude, & - model_state_variables, & - debug, & - grid_oi - -contains - -!------------------------------------------------------------------ -! -! Called to do one time initialization of the model. As examples, -! might define information about the model size or model timestep. -! In models that require pre-computed static data, for instance -! spherical harmonic weights, these would also be computed here. -! Can be a NULL INTERFACE for the simplest models. - -subroutine static_init_model() - - real(r8) :: x_loc - integer :: i, dom_id,iunit,io,ss,dd -!integer :: iunit, io - -if ( module_initialized ) return ! only need to do this once. - -! Print module information to log file and stdout. -call register_module(source, revision, revdate) - -module_initialized = .true. - -! This is where you would read a namelist, for example. -call find_namelist_in_file("input.nml", "model_nml", iunit) -read(iunit, nml = model_nml, iostat = io) -call check_namelist_read(iunit, io, "model_nml") - -call error_handler(E_MSG,'static_init_model','model_nml values are',' ',' ',' ') -if (do_nml_file()) write(nmlfileunit, nml=model_nml) -if (do_nml_term()) write( * , nml=model_nml) - -call set_calendar_type('Gregorian') - -model_timestep = set_model_time_step() - -call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400) - -write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds' -call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate) - -call get_horiz_grid_dims(Nx) -call get_ncat_dim(Ncat) - -call verify_state_variables(model_state_variables, nfields, variable_table, & - state_kinds_list, update_var_list) - -allocate(TLAT(Nx), TLON(Nx)) - -call read_horiz_grid(Nx, TLAT, TLON) - -if (do_output()) write(logfileunit, *) 'Using grid : Nx, Ncat = ', & - Nx, Ncat -if (do_output()) write( * , *) 'Using grid : Nx, Ncat = ', & - Nx, Ncat - -domain_id = add_domain('cice.r.nc', nfields, & - var_names = variable_table(1:nfields, VAR_NAME_INDEX), & - kind_list = state_kinds_list(1:nfields), & - update_list = update_var_list(1:nfields)) - -if (debug > 2) call state_structure_info(domain_id) - -model_size = get_domain_size(domain_id) -if (do_output()) write(*,*) 'model_size = ', model_size - - -end subroutine static_init_model -!------------------------------------------------------------------ -! Returns a model state vector, x, that is some sort of appropriate -! initial condition for starting up a long integration of the model. -! At present, this is only used if the namelist parameter -! start_from_restart is set to .false. in the program perfect_model_obs. -! If this option is not to be used in perfect_model_obs, or if no -! synthetic data experiments using perfect_model_obs are planned, -! this can be a NULL INTERFACE. - -!subroutine init_conditions(x) -! -!real(r8), intent(out) :: x(:) -! -!x = MISSING_R8 -! -!end subroutine init_conditions - - - -!------------------------------------------------------------------ -! Does a single timestep advance of the model. The input value of -! the vector x is the starting condition and x is updated to reflect -! the changed state after a timestep. The time argument is intent -! in and is used for models that need to know the date/time to -! compute a timestep, for instance for radiation computations. -! This interface is only called if the namelist parameter -! async is set to 0 in perfect_model_obs of filter or if the -! program integrate_model is to be used to advance the model -! state as a separate executable. If one of these options -! is not going to be used (the model will only be advanced as -! a separate model-specific executable), this can be a -! NULL INTERFACE. - -!subroutine adv_1step(x, time) -! -!real(r8), intent(inout) :: x(:) -!type(time_type), intent(in) :: time -! -!end subroutine adv_1step - - - -!------------------------------------------------------------------ -! Returns the number of items in the state vector as an integer. -! This interface is required for all applications. - -function get_model_size() - -integer(i8) :: get_model_size - -get_model_size = model_size - -end function get_model_size - - - -!------------------------------------------------------------------ -! Companion interface to init_conditions. Returns a time that is somehow -! appropriate for starting up a long integration of the model. -! At present, this is only used if the namelist parameter -! start_from_restart is set to .false. in the program perfect_model_obs. -! If this option is not to be used in perfect_model_obs, or if no -! synthetic data experiments using perfect_model_obs are planned, -! this can be a NULL INTERFACE. - -!subroutine init_time(time) -! -!type(time_type), intent(out) :: time -! -!! for now, just set to 0 -!time = set_time(0,0) -! -!end subroutine init_time - -!------------------------------------------------------------------ -! Given a state handle, a location, and a model state variable type, -! interpolates the state variable fields to that location and returns -! the values in expected_obs. The istatus variables should be returned as -! 0 unless there is some problem in computing the interpolation in -! which case an alternate value should be returned. The itype variable -! is a model specific integer that specifies the kind of field (for -! instance temperature, zonal wind component, etc.). In low order -! models that have no notion of types of variables this argument can -! be ignored. For applications in which only perfect model experiments -! with identity observations (i.e. only the value of a particular -! state variable is observed), this can be a NULL INTERFACE. - -subroutine model_interpolate(state_handle, ens_size, location, obs_type, expected_obs, istatus, thick_flag) - - -type(ensemble_type), intent(in) :: state_handle -integer, intent(in) :: ens_size -type(location_type), intent(in) :: location -integer, intent(in) :: obs_type -real(r8), intent(out) :: expected_obs(ens_size) !< array of interpolated values -integer, intent(out) :: istatus(ens_size) -logical,optional, intent(inout) :: thick_flag - -!local vars -real(r8) :: loc_array(3), llon, llat -integer(i8) :: base_offset -integer :: cat_index, cat_signal, icat, cat_signal_interm -real(r8) :: expected_aggr_conc(ens_size) -integer :: set_obstype -integer :: var_table_index - -!Fei---need aicen*fyn to calculate the aggregate FY concentration------------ -real(r8) :: expected_conc(ens_size) -real(r8) :: expected_fy(ens_size) -real(r8) :: expected_tsfc(ens_size) -real(r8) :: temp(ens_size) -real(r8) :: temp1(ens_size) - -if ( .not. module_initialized ) call static_init_model - -expected_obs(:) = MISSING_R8 ! the DART bad value flag -istatus(:) = 99 - -loc_array = get_location(location) -llon = loc_array(1) -llat = loc_array(2) -cat_index = int(loc_array(3)) - -if (obs_type == QTY_SEAICE_CATEGORY) then - if (cat_index <= Ncat) then - istatus = 0 - expected_obs = cat_index - RETURN - endif -endif -if (debug > 1) then - print *, 'requesting interpolation of ', obs_type, ' at ', llon, llat, cat_index -endif - -SELECT CASE (obs_type) - CASE (QTY_SEAICE_AGREG_THICKNESS ) ! these kinds require aggregating 3D vars to make a 2D var - if (any(variable_table(:,1)=='hi')) then - cat_signal = 1 !was 1 ! for extra special procedure to aggregate - !base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_AGREG_THICKNESS)) - thick_flag = .true. - base_offset = cat_index - set_obstype = obs_type - !call find_var_type('hi',var_index) - else - set_obstype = QTY_SEAICE_VOLUME - cat_signal = 1 ! for extra special procedure to aggregate - !base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_VOLUME)) - base_offset = cat_index - !call find_var_type('vicen',var_index) - endif - CASE (QTY_SEAICE_AGREG_SNOWDEPTH ) ! these kinds require aggregating 3D vars to make a 2D var - if (any(variable_table(:,1)=='hs')) then - cat_signal = 1 !was 1 ! for extra special procedure to aggregate - !base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_AGREG_SNOWDEPTH)) - base_offset = cat_index - thick_flag = .true. - set_obstype = obs_type - !call find_var_type('hs',var_index) - else - set_obstype = QTY_SEAICE_SNOWVOLUME - cat_signal = 1 ! for extra special procedure to aggregate - !base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SNOWVOLUME)) - base_offset = cat_index - !call find_var_type('vsnon',var_index) - endif - CASE (QTY_SEAICE_AGREG_CONCENTR ) ! these kinds require aggregating a 3D var to make a 2D var - cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - set_obstype = obs_type - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_CONCENTR)) - CASE (QTY_SEAICE_AGREG_VOLUME ) ! these kinds require aggregating a 3D var to make a 2D var - cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - set_obstype = obs_type - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_VOLUME)) - CASE (QTY_SEAICE_AGREG_SNOWVOLUME ) ! these kinds require aggregating a 3D var to make a 2D var - cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - set_obstype = obs_type - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SNOWVOLUME)) - CASE (QTY_SEAICE_AGREG_SURFACETEMP) ! FEI need aicen to average the temp, have not considered open water temp yet - if (any(variable_table(:,1)=='Tsfc')) then - cat_signal = 1 - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_AGREG_SURFACETEMP)) - thick_flag = .true. - set_obstype = obs_type - else - cat_signal = -3 - set_obstype = QTY_SEAICE_SURFACETEMP - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SURFACETEMP)) - endif - CASE (QTY_SOM_TEMPERATURE) ! these kinds are 1d variables - cat_signal = 1 - set_obstype = obs_type - !base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SOM_TEMPERATURE)) - base_offset = cat_index - CASE (QTY_SEAICE_CONCENTR , & ! these kinds have an additional dim for category - QTY_SEAICE_FY , & - QTY_SEAICE_VOLUME , & - QTY_SEAICE_SNOWVOLUME , & - QTY_SEAICE_SURFACETEMP , & - QTY_SEAICE_FIRSTYEARAREA , & - QTY_SEAICE_ICEAGE , & - QTY_SEAICE_LEVELAREA , & - QTY_SEAICE_LEVELVOLUME , & - QTY_SEAICE_MELTPONDAREA , & - QTY_SEAICE_MELTPONDDEPTH , & - QTY_SEAICE_MELTPONDLID , & - QTY_SEAICE_MELTPONDSNOW , & - QTY_SEAICE_SALINITY001 , & - QTY_SEAICE_SALINITY002 , & - QTY_SEAICE_SALINITY003 , & - QTY_SEAICE_SALINITY004 , & - QTY_SEAICE_SALINITY005 , & - QTY_SEAICE_SALINITY006 , & - QTY_SEAICE_SALINITY007 , & - QTY_SEAICE_SALINITY008 , & - QTY_SEAICE_ICEENTHALPY001 , & - QTY_SEAICE_ICEENTHALPY002 , & - QTY_SEAICE_ICEENTHALPY003 , & - QTY_SEAICE_ICEENTHALPY004 , & - QTY_SEAICE_ICEENTHALPY005 , & - QTY_SEAICE_ICEENTHALPY006 , & - QTY_SEAICE_ICEENTHALPY007 , & - QTY_SEAICE_ICEENTHALPY008 , & - QTY_SEAICE_SNOWENTHALPY001, & - QTY_SEAICE_SNOWENTHALPY002, & - QTY_SEAICE_SNOWENTHALPY003 ) - ! move pointer to the particular category - ! then treat as 2d field in lon_lat_interp - - base_offset = get_index_start(domain_id, get_varid_from_kind(obs_type)) - base_offset = base_offset + (cat_index-1)! * Nx - base_offset = cat_index - set_obstype = obs_type - cat_signal = 1 ! now same as boring 2d field - CASE DEFAULT - ! Not a legal type for interpolation, return istatus error - istatus = 15 - return -END SELECT - -if (cat_signal == -2) then - temp = 0.0_r8 - temp1= 0.0_r8 - do icat = 1,Ncat - !reads in aicen - cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_CONCENTR)) - base_offset = base_offset + (icat-1) * Nx - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal_interm, expected_conc, istatus) - !reads in fyn - cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_FY)) - base_offset = base_offset + (icat-1) * Nx - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal_interm, expected_fy, istatus) - temp = temp + expected_conc * expected_fy !sum(aicen*fyn) = FY % over ice - temp1= temp1+ expected_conc !sum(aicen) = aice - - if (any(expected_conc<0.0) .or. any(expected_conc>1.0))then - print*,'obstype FY expected sicn:',expected_conc - print*,'FY sicn lat lon:',llat,llon - endif - if (any(expected_fy>1.0) .or. any(expected_fy<0.0)) then - print*,'obstype FY expected fyn:',expected_fy,llat,llon - print*,'FY fyn lat lon:',llat,llon - endif - - end do - expected_obs = temp/max(temp1,1.0e-8) !sum(aicen*fyn)/aice = FY % in the gridcell -else if (cat_signal == -3 ) then - temp = 0.0_r8 - temp1= 0.0_r8 - do icat = 1,Ncat - !reads in aicen - cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_CONCENTR)) - base_offset = base_offset + (icat-1) * Nx - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal_interm, expected_conc, istatus) - !reads in Tsfcn - cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_SURFACETEMP)) - base_offset = base_offset + (icat-1) * Nx - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal_interm, expected_tsfc, istatus) - if (any(expected_conc<0.0) .or. any(expected_conc>1.0))then - print*,'obstype TSFC expected sicn:',expected_conc - print*,'TSFC sicn lat lon:',llat,llon - endif - if (any(expected_tsfc>50.0) .or. any(expected_tsfc<-100.0)) then - print*,'obstype TSFC expected tsfcn:',expected_tsfc - print*,'TSFC tsfcn lat lon:',llat,llon - endif - temp = temp + expected_conc * expected_tsfc !sum(aicen*Tsfcn) - temp1= temp1+ expected_conc !sum(aicen) = aice - end do - expected_obs = temp/max(temp1,1.0e-8) !sum(aicen*Tsfcn)/aice = Tsfc ;averaged temperature over sea-ice covered portion - if (any(expected_obs>50.0) .or. any(expected_obs<-100.0)) then - print*,'obstype TSFC expected obs:',expected_obs - print*,'TSFC tsfc lat lon:' ,llat,llon - print*,'temp:',temp - print*,'temp1:',temp1 - endif -else - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal, expected_obs, istatus) - - if (any(expected_obs<0.0))then - print*,'obstype SIC expected concs:',expected_obs - print*,'SIC sic negative lat lon:',llat,llon - endif - if (any(expected_obs>1.0))then - print*,'obstype SIC expected concs:',expected_obs - print*,'SIC sic positive lat lon:',llat,llon - endif -endif - -if (cat_signal == -1) then - ! we need to know the aggregate sea ice concentration for these special cases - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_CONCENTR)) - base_offset = base_offset + (cat_index-1) - print*,'CHECK CHECK CHECK' - call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, set_obstype, cat_signal, expected_aggr_conc, istatus) - expected_obs = expected_obs/max(expected_aggr_conc,1.0e-8) ! hope this is allowed so we never divide by zero - - if (any(expected_aggr_conc<0.0) .or. any(expected_aggr_conc>1.0))then - print*,'obstype SIT expected conc:',expected_aggr_conc - print*,'SIT sic lat lon:',llat,llon - endif - -endif - -if (debug > 1) print *, 'interp val, istatus = ', expected_obs, istatus, size(expected_obs) - -! This should be the result of the interpolation of a -! given kind (itype) of variable at the given location. - -! The return code for successful return should be 0. -! Any positive number is an error. -! Negative values are reserved for use by the DART framework. -! Using distinct positive values for different types of errors can be -! useful in diagnosing problems. - -end subroutine model_interpolate -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine lon_lat_interpolate(state_handle, ens_size, offset, lon, lat, var_type, cat_signal, expected_obs, istatus) -type(ensemble_type), intent(in) :: state_handle -integer, intent(in) :: ens_size -integer(i8), intent(in) :: offset -real(r8), intent(in) :: lon, lat -integer, intent(in) :: var_type -integer, intent(in) :: cat_signal -real(r8), intent(out) :: expected_obs(ens_size) -integer, intent(out) :: istatus(ens_size) - -integer :: lat_bot, lat_top, lon_bot, lon_top, num_inds, start_ind -integer :: x_ind, y_ind -real(r8) :: x_corners(4), y_corners(4) -real(r8) :: p(4,ens_size), xbot(ens_size), xtop(ens_size) -real(r8) :: work_expected_obs(ens_size) -real(r8) :: lon_fract, lat_fract -logical :: masked -integer :: quad_status -integer :: e, iterations, Niterations -integer :: next_offset -integer(i8) :: state_index -if ( .not. module_initialized ) call static_init_model - -istatus = 0 -print*,'VAR TYPE',var_type -if (var_type == 14) then - e = 1 -else if (var_type == 15) then - e = 2 -else if (var_type == 16) then - e = 3 -endif -if ( cat_signal < 1 ) then - Niterations = Ncat ! only iterate if aggregating over all types -else - Niterations = 1 ! no need to iterate -endif -work_expected_obs = 0.0_r8 -expected_obs = 0.0_r8 -do iterations = 1, Niterations - - ! FIXME: this should use the state structure routine 'get_dart_vector_index' - ! to get the start of the next category layer. this code assumes it knows - ! exactly how the state vector is laid out (reasonable, but might not be true - ! in future versions of dart.) - !next_offset = offset + (iterations-1)*Nx - !print*,'offset',offset - state_index = get_dart_vector_index(grid_oi,int(offset,i4),1, domain_id, e) - work_expected_obs = get_state(state_index,state_handle) - !if(masked) then - ! istatus = 3 - ! return - !endif - expected_obs = expected_obs+work_expected_obs -enddo -end subroutine lon_lat_interpolate -!------------------------------------------------------------------ -! Returns the smallest increment in time that the model is capable -! of advancing the state in a given implementation, or the shortest -! time you want the model to advance between assimilations. -! This interface is required for all applications. - -function shortest_time_between_assimilations() - -type(time_type) :: shortest_time_between_assimilations - -if ( .not. module_initialized ) call static_init_model - -shortest_time_between_assimilations = model_timestep - -end function shortest_time_between_assimilations -!------------------------------------------------------------------ -! Given an integer index into the state vector structure, returns the -! associated location. A second intent(out) optional argument kind -! can be returned if the model has more than one type of field (for -! instance temperature and zonal wind component). This interface is -! required for all filter applications as it is required for computing -! the distance between observations and state variables. - -subroutine get_state_meta_data(index_in, location, var_type) - -integer(i8), intent(in) :: index_in -type(location_type), intent(out) :: location -integer, intent(out), optional :: var_type - -real(r8) :: lat, lon, rcat -integer :: ni_index, hold_index, cat_index, local_var, var_id - -! these should be set to the actual location and state quantity -if ( .not. module_initialized ) call static_init_model - -call get_model_variable_indices(index_in, ni_index, cat_index, hold_index, var_id=var_id) -call get_state_kind(var_id, local_var) - -lon = TLON(ni_index) -lat = TLAT(ni_index) - -if (debug > 5) print *, 'lon, lat, cat_index = ', lon, lat, cat_index -rcat = cat_index*1.0_r8 -location = set_location(lon, lat, rcat, VERTISLEVEL) - -if (present(var_type)) then - var_type = local_var -endif - -end subroutine get_state_meta_data - -subroutine get_state_kind(var_ind, var_type) - integer, intent(in) :: var_ind - integer, intent(out) :: var_type - -! Given an integer index into the state vector structure, returns the kind, -! and both the starting offset for this kind, as well as the offset into -! the block of this kind. - -if ( .not. module_initialized ) call static_init_model - -var_type = state_kinds_list(var_ind) - -end subroutine get_state_kind - - -!------------------------------------------------------------------ -! Does any shutdown and clean-up needed for model. Can be a NULL -! INTERFACE if the model has no need to clean up storage, etc. - -subroutine end_model() - -deallocate(TLAT,TLON) - -end subroutine end_model - - -!------------------------------------------------------------------ -! write any additional attributes to the output and diagnostic files - -subroutine nc_write_model_atts(ncid, domain_id) - -integer, intent(in) :: ncid ! netCDF file identifier -integer, intent(in) :: domain_id -integer :: NGridDimID - -integer, parameter :: MAXLINELEN = 128 -character(len=8), parameter :: cice_namelist_file = 'cice_in' -character(len=MAXLINELEN), allocatable, dimension(:) :: textblock -integer :: LineLenDimID, nlinesDimID, nmlVarID -integer :: nlines, linelen,status -logical :: has_cice_namelist - -character(len=256) :: filename - -integer :: NlonDimID, NlatDimID -integer :: tlonVarID, tlatVarID - -if ( .not. module_initialized ) call static_init_model - -! put file into define mode. - -write(filename,*) 'ncid', ncid - -call nc_begin_define_mode(ncid) - -call nc_add_global_creation_time(ncid) - -call nc_add_global_creation_time(ncid) - -call nc_add_global_attribute(ncid, "model_source", source ) -call nc_add_global_attribute(ncid, "model_revision", revision ) -call nc_add_global_attribute(ncid, "model_revdate", revdate ) -call nc_add_global_attribute(ncid, "model", "CICE-SCM") - -call nc_check(nf90_def_dim(ncid, name='ni', & - len = Nx, dimid = NGridDimID),'nc_write_model_atts', 'ni def_dim '//trim(filename)) - -call nc_check(nf90_def_var(ncid,name='TLON', xtype=nf90_real, & - dimids=(/ NGridDimID /), varid=tlonVarID),& - 'nc_write_model_atts', 'TLON def_var '//trim(filename)) -call nc_check(nf90_def_var(ncid,name='TLAT', xtype=nf90_real, & - dimids=(/ NGridDimID /), varid=tlatVarID),& - 'nc_write_model_atts', 'TLAT def_var '//trim(filename)) - -call nc_end_define_mode(ncid) - -call nc_check(nf90_put_var(ncid, tlonVarID, TLON ), & - 'nc_write_model_atts', 'TLON put_var '//trim(filename)) -call nc_check(nf90_put_var(ncid, tlatVarID, TLAT ), & - 'nc_write_model_atts', 'TLAT put_var '//trim(filename)) - -! Flush the buffer and leave netCDF file open -call nc_synchronize_file(ncid) - -end subroutine nc_write_model_atts -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -function get_varid_from_kind(dart_kind) - -integer, intent(in) :: dart_kind -integer :: get_varid_from_kind - -! given a kind, return what variable number it is - -integer :: i - -do i = 1, get_num_variables(domain_id) - if (dart_kind == state_kinds_list(i)) then - get_varid_from_kind = i - return - endif -end do - -if (debug > 1) then - write(string1, *) 'Kind ', dart_kind, ' not found in state vector' - write(string2, *) 'AKA ', get_name_for_quantity(dart_kind), ' not found in state vector' - call error_handler(E_MSG,'get_varid_from_kind', string1, & - source, revision, revdate, text2=string2) -endif - -get_varid_from_kind = -1 - -end function get_varid_from_kind -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine verify_state_variables( state_variables, ngood, table, kind_list, update_var ) - -character(len=*), intent(inout) :: state_variables(:) -integer, intent(out) :: ngood -character(len=*), intent(out) :: table(:,:) -integer, intent(out) :: kind_list(:) ! kind number -logical, optional, intent(out) :: update_var(:) ! logical update - -integer :: nrows, i -character(len=NF90_MAX_NAME) :: varname, dartstr, update - -if ( .not. module_initialized ) call static_init_model - -nrows = size(table,1) - -ngood = 0 - -!>@todo deprecate. Remove a hidden 'default' set of variables. -!>@ The default is provided in the input namelist. - -if ( state_variables(1) == ' ' ) then ! no model_state_variables namelist provided - call use_default_state_variables( state_variables ) - string1 = 'model_nml:model_state_variables not specified using default variables' - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) -endif - -MyLoop : do i = 1, nrows - - varname = trim(state_variables(3*i -2)) - dartstr = trim(state_variables(3*i -1)) - update = trim(state_variables(3*i )) - - call to_upper(update) - - table(i,1) = trim(varname) - table(i,2) = trim(dartstr) - table(i,3) = trim(update) - - if ( table(i,1) == ' ' .and. table(i,2) == ' ' .and. table(i,3) == ' ') exit MyLoop - - if ( table(i,1) == ' ' .or. table(i,2) == ' ' .or. table(i,3) == ' ' ) then - string1 = 'model_nml:model_state_variables not fully specified' - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure DART kind is valid - - kind_list(i) = get_index_for_quantity(dartstr) - if( kind_list(i) < 0 ) then - write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure the update variable has a valid name - - if ( present(update_var) )then - SELECT CASE (update) - CASE ('UPDATE') - update_var(i) = .true. - CASE ('NO_COPY_BACK') - update_var(i) = .false. - CASE DEFAULT - write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' - write(string2,'(6A)') 'you provided : ', trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate, text2=string2) - END SELECT - endif - - ! Record the contents of the DART state vector - - if (do_output()) then - write(string1,'(A,I2,6A)') 'variable ',i,' is ',trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) - endif - - ngood = ngood + 1 -enddo MyLoop - -end subroutine verify_state_variables -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine use_default_state_variables( state_variables ) - -character(len=*), intent(inout) :: state_variables(:) - -! strings must all be the same length for the gnu compiler -state_variables( 1:5*num_state_table_columns ) = & - (/ 'CONCENTRATION ', 'QTY_SEAICE_CONCENTR ', 'UPDATE ', & - 'ICEVOLUME ', 'QTY_SEAICE_VOLUME ', 'UPDATE ', & - 'SNOWVOLUME ', 'QTY_SEAICE_SNOWVOLUME ', 'UPDATE ', & - 'UICE ', 'QTY_U_SEAICE_COMPONENT ', 'UPDATE ', & - 'VICE ', 'QTY_V_SEAICE_COMPONENT ', 'UPDATE '/) - -end subroutine use_default_state_variables -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine get_close_state(filt_gc, base_loc, base_type, locs, loc_qtys, loc_indx, & - num_close, close_indices, distances, state_handle) - -type(get_close_type), intent(in) :: filt_gc -type(location_type), intent(inout) :: base_loc -integer, intent(in) :: base_type -type(location_type), intent(inout) :: locs(:) -integer, intent(in) :: loc_qtys(:) -integer(i8), intent(in) :: loc_indx(:) -integer, intent(out) :: num_close -integer, intent(out) :: close_indices(:) -real(r8), intent(out), optional :: distances(:) -type(ensemble_type), intent(in), optional :: state_handle - -! Given a DART location (referred to as "base") and a set of candidate -! locations & kinds (locs, loc_qtys/indx), returns the subset close to the -! "base", their indices, and their distances to the "base" ... - -integer :: t_ind, k - -! Initialize variables to missing status - -num_close = 0 -close_indices = -99 -if (present(distances)) distances(:) = 1.0e9 !something big and positive (far away) - -! Get all the potentially close obs but no dist (optional argument dist(:) -! is not present) This way, we are decreasing the number of distance -! computations that will follow. This is a horizontal-distance operation and -! we don't need to have the relevant vertical coordinate information yet -! (for obs). -call loc_get_close_state(filt_gc, base_loc, base_type, locs, loc_qtys, loc_indx, & - num_close, close_indices) - -! Loop over potentially close subset of obs priors or state variables -if (present(distances)) then - do k = 1, num_close - - t_ind = close_indices(k) - - ! if dry land, leave original 1e9 value. otherwise, compute real dist. - distances(k) = get_dist(base_loc, locs(t_ind), & - base_type, loc_qtys(t_ind)) - enddo -endif - -end subroutine get_close_state -!!!!!!!!!!!!!!!! -function read_model_time(filename) - -character(len=256) :: filename -type(time_type) :: read_model_time - -integer :: ncid !< netcdf file id -integer :: nyr , & ! year number, in cice restart - month , & ! month number, 1 to 12, in cice restart - mday , & ! day of the month, in cice restart - sec ! elapsed seconds into date, in cice restart -integer :: hour , & ! hour of the day, needed for dart set_date - minute , & ! minute of the hour, needed for dart set_date - secthismin - -if ( .not. module_initialized ) call static_init_model - -if ( .not. file_exist(filename) ) then - write(string1,*) 'cannot open file ', trim(filename),' for reading.' - call error_handler(E_ERR,'read_model_time',string1,source,revision,revdate) -endif - -call nc_check( nf90_open(trim(filename), NF90_NOWRITE, ncid), & - 'read_model_time', 'open '//trim(filename)) -call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'nyr' , nyr), & - 'read_model_time', 'get_att nyr') -call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'month' , month), & - 'read_model_time', 'get_att month') -call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'mday' , mday), & - 'read_model_time', 'get_att mday') -call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'sec', sec), & - 'read_model_time', 'get_att sec') - -! FIXME: we don't allow a real year of 0 - add one for now, but -! THIS MUST BE FIXED IN ANOTHER WAY! -if (nyr == 0) then - call error_handler(E_MSG, 'read_model_time', & - 'WARNING!!! year 0 not supported; setting to year 1') - nyr = 1 -endif - -hour = int(sec/3600) -minute = int((sec-hour*3600)/60) -secthismin = int(sec-hour*3600-minute*60) - -read_model_time = set_date(nyr, month, mday, hour, minute, secthismin) -end function read_model_time -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine write_model_time(ncid, model_time, adv_to_time) - -integer, intent(in) :: ncid -type(time_type), intent(in) :: model_time -type(time_type), intent(in), optional :: adv_to_time - -character(len=16), parameter :: routine = 'write_model_time' - -integer :: io, varid, iyear, imonth, iday, ihour, imin, isec -integer :: seconds - -if ( .not. module_initialized ) call static_init_model - -if (present(adv_to_time)) then - call get_date(adv_to_time, iyear, imonth, iday, ihour, imin, isec) - write(string1,*)'CICE/DART not configured to advance CICE.' - write(string2,*)'called with optional advance_to_time of' - write(string3,'(i4.4,5(1x,i2.2))')iyear,imonth,iday,ihour,imin, isec - call error_handler(E_ERR, routine, string1, & - source, revision, revdate, text2=string2,text3=string3) -endif - -call get_date(model_time, iyear, imonth, iday, ihour, imin, isec) - -seconds = (ihour*60 + imin)*60 + isec - -call nc_begin_define_mode(ncid) -call nc_add_global_attribute(ncid, 'nyr' , iyear) -call nc_add_global_attribute(ncid, 'month' , imonth) -call nc_add_global_attribute(ncid, 'mday' , iday) -call nc_add_global_attribute(ncid, 'sec' , seconds) -call nc_end_define_mode(ncid) - -end subroutine write_model_time -!----------------------------------------------------------------- -! Check which surface temperature state variable is in restart -subroutine check_sfctemp_var(flag) -logical, intent(inout) :: flag - -if (any(variable_table(:,1)=='Tsfc')) then - flag = .true. -else - flag = .false. -endif -end subroutine check_sfctemp_var -!----------------------------------------------------------------- -! Find state variable index -subroutine find_var_type(varname,var_index) -character(len=16), intent(in) :: varname -integer, intent(inout) :: var_index - -integer :: i - -do i=1,size(variable_table(:,1)) - if (trim(varname) == variable_table(i,1)) then - var_index = i - return - endif -enddo -write(string1,*)'Could not find index of state variable' -call error_handler(E_ERR, 'find_var_type', string1, & - source, revision, revdate) -end subroutine find_var_type -!=================================================================== -! End of model_mod -!=================================================================== -end module model_mod - -! -! $URL$ -! $Id$ -! $Revision$ -! $Date$ diff --git a/models/cice-scm2/readme.rst b/models/cice-scm2/readme.rst deleted file mode 100644 index 1b867a5daf..0000000000 --- a/models/cice-scm2/readme.rst +++ /dev/null @@ -1,5 +0,0 @@ -cice-scm2 -============== - -.. attention:: - Add your model documentation here. diff --git a/models/cice-scm2/work/algorithm_info_mod.f90 b/models/cice-scm2/work/algorithm_info_mod.f90 deleted file mode 100644 index 19bb14e1f9..0000000000 --- a/models/cice-scm2/work/algorithm_info_mod.f90 +++ /dev/null @@ -1,215 +0,0 @@ -! DART software - Copyright UCAR. This open source software is provided -! by UCAR, "as is", without charge, subject to all terms of use at -! http://www.image.ucar.edu/DAReS/DART/DART_download - -module algorithm_info_mod - -use types_mod, only : r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_SEAICE_VOLUME, QTY_SEAICE_CONCENTR, QTY_SEAICE_SNOWVOLUME, & - QTY_SEAICE_AGREG_THICKNESS, QTY_SEAICE_AGREG_CONCENTR, QTY_SEAICE_AGREG_FREEBOARD -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata - -implicit none -private - -integer, parameter :: NORMAL_PRIOR = 1 -integer, parameter :: BOUNDED_NORMAL_RH_PRIOR = 2 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - NORMAL_PRIOR, BOUNDED_NORMAL_RH_PRIOR - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, bounded, bounds) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded(2) -real(r8), intent(out) :: bounds(2) - -integer :: obs_type, obs_kind - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -obs_kind = get_quantity_for_type_of_obs(obs_type) - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -if(obs_kind == QTY_SEAICE_AGREG_CONCENTR) then - bounded(1) = .true.; bounded(2) = .true. - bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 -elseif(obs_kind == QTY_SEAICE_AGREG_THICKNESS) then - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; -elseif(obs_kind == QTY_SEAICE_AGREG_FREEBOARD) then - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; -else - bounded = .false. -endif - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded, bounds) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded(2) -real(r8), intent(out) :: bounds(2) - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_PRIOR -! which does nothing or BOUNDED_NORMAL_RH_PRIOR. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded(1) = .true.; bounded(2) = .true. -! bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -if(is_inflation) then - ! Case for inflation transformation - if(kind == QTY_SEAICE_CONCENTR) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .true. - bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 - elseif(kind == QTY_SEAICE_VOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - elseif(kind == QTY_SEAICE_SNOWVOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - else - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded = .false. - endif -elseif(is_state) then - ! Case for state variable priors - if(kind == QTY_SEAICE_CONCENTR) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .true. - bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 - elseif(kind == QTY_SEAICE_VOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - elseif(kind == QTY_SEAICE_SNOWVOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - else - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded = .false. - endif -else - ! This case is for observation (extended state) priors - if(kind == QTY_SEAICE_CONCENTR) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .true. - bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 - elseif(kind == QTY_SEAICE_VOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - elseif(kind == QTY_SEAICE_SNOWVOLUME) then - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; - else - dist_type = BOUNDED_NORMAL_RH_PRIOR - bounded = .false. - endif -endif - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded, bounds) - -integer, intent(in) :: obs_kind -integer, intent(out) :: filter_kind -logical, intent(out) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(out) :: sort_obs_inc -logical, intent(out) :: spread_restoration -logical, intent(out) :: bounded(2) -real(r8), intent(out) :: bounds(2) - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_SEAICE_AGREG_CONCENTR) then - filter_kind = 101 - bounded(1) = .true.; bounded(2) = .true. - bounds(1) = 0.0_r8; bounds(2) = 1.0_r8 -elseif(obs_kind == QTY_SEAICE_AGREG_THICKNESS) then - filter_kind = 101 - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; -elseif(obs_kind == QTY_SEAICE_AGREG_FREEBOARD) then - filter_kind = 101 - bounded(1) = .true.; bounded(2) = .false. - bounds(1) = 0.0_r8; -else - filter_kind = 101 - bounded = .false. -endif - -! HK you are overwritting filter kind in the if statement with this: -filter_kind = 101 - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options on old RHF implementation -! rectangular_quadrature = .true. -! gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/models/cice-scm2/work/input.nml b/models/cice-scm2/work/input.nml deleted file mode 100644 index 56706e4141..0000000000 --- a/models/cice-scm2/work/input.nml +++ /dev/null @@ -1,220 +0,0 @@ -&perfect_model_obs_nml - read_input_state_from_file = .true., - single_file_in = .false. - input_state_files = "input_file.nc" - - write_output_state_to_file = .false., - single_file_out = .true. - output_state_files = "perfect_output.nc" - output_interval = 1, - - async = 0, - adv_ens_command = "./advance_model.csh", - - obs_seq_in_file_name = "obs_seq.in", - obs_seq_out_file_name = "obs_seq.out", - init_time_days = 0, - init_time_seconds = 0, - first_obs_days = -1, - first_obs_seconds = -1, - last_obs_days = -1, - last_obs_seconds = -1, - - trace_execution = .false., - output_timestamps = .false., - print_every_nth_obs = -1, - output_forward_op_errors = .false., - silence = .false., - / - -&filter_nml - single_file_in = .true., - input_state_files = '' - input_state_file_list = 'filter_input_list.txt' - - stages_to_write = 'input', 'preassim', 'analysis', 'output' - - single_file_out = .true., - output_state_files = '' - output_state_file_list = 'filter_output_list.txt' - output_interval = 1, - output_members = .true. - num_output_state_members = 0, - output_mean = .true. - output_sd = .true. - write_all_stages_at_end = .false. - - ens_size = 29, - num_groups = 1, - perturb_from_single_instance = .false., - perturbation_amplitude = 0.2, - distributed_state = .true. - - async = 0, - adv_ens_command = "./advance_model.csh", - - obs_sequence_in_name = "obs_seq.out", - obs_sequence_out_name = "obs_seq.final", - num_output_obs_members = 20, - init_time_days = 0, - init_time_seconds = 0, - first_obs_days = -1, - first_obs_seconds = -1, - last_obs_days = -1, - last_obs_seconds = -1, - - inf_flavor = 0, 0, - inf_initial_from_restart = .false., .false., - inf_sd_initial_from_restart = .false., .false., - inf_deterministic = .true., .true., - inf_initial = 1.0, 1.0, - inf_lower_bound = 1.0, 1.0, - inf_upper_bound = 100.0, 1000000.0, - inf_damping = 1.0, 1.0, - inf_sd_initial = 0.0, 0.0, - inf_sd_lower_bound = 0.0, 0.0, - inf_sd_max_change = 1.05, 1.05, - - trace_execution = .false., - output_timestamps = .false., - output_forward_op_errors = .false., - silence = .false., - / - -&smoother_nml - num_lags = 0, - start_from_restart = .false., - output_restart = .false., - restart_in_file_name = 'smoother_ics', - restart_out_file_name = 'smoother_restart' - / - -&ensemble_manager_nml - / - -&assim_tools_nml - filter_kind = 1, - cutoff = 1000000.0 - sort_obs_inc = .false., - spread_restoration = .false., - sampling_error_correction = .false., - adaptive_localization_threshold = -1, - distribute_mean = .false. - output_localization_diagnostics = .false., - localization_diagnostics_file = 'localization_diagnostics', - print_every_nth_obs = 0 - / - -&cov_cutoff_nml - select_localization = 1 - / - -®_factor_nml - select_regression = 1, - input_reg_file = "time_mean_reg", - save_reg_diagnostics = .false., - reg_diagnostics_file = "reg_diagnostics" - / - -&obs_sequence_nml - write_binary_obs_sequence = .false. - / - -&obs_kind_nml - assimilate_these_obs_types = 'SAT_SEAICE_AGREG_THICKNESS' - evaluate_these_obs_types = '' - / - -&model_nml - assimilation_period_days = 1 - assimilation_period_seconds = 0 - model_perturbation_amplitude = 2e-05 - debug = 100 - model_state_variables = 'aicen', 'QTY_SEAICE_CONCENTR', 'UPDATE', 'vicen', - 'QTY_SEAICE_VOLUME', 'UPDATE', 'vsnon', 'QTY_SEAICE_SNOWVOLUME', - 'UPDATE' -/ - -&dart_to_cice_nml - dart_to_cice_input_file = 'restart_state.nc' - original_cice_input_file = 'dart_restart.nc' - previous_cice_input_file = 'pre_restart.nc' - balance_method = 'simple_squeeze' - r_snw_name = 'r_snw_vary' - gridpt_oi = 3 -/ - -&utilities_nml - TERMLEVEL = 1, - module_details = .false., - logfilename = 'dart_log.out', - nmlfilename = 'dart_log.nml', - write_nml = 'none' - / - -&preprocess_nml - input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' - output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' - input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' - output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' - obs_type_files = '../../../observations/forward_operators/obs_def_cice_mod.f90' - quantity_files = '../../../assimilation_code/modules/observations/seaice_quantities_mod.f90', - '../../../assimilation_code/modules/observations/ocean_quantities_mod.f90' - / - -&obs_sequence_tool_nml - filename_seq = 'obs_seq.one', 'obs_seq.two', - filename_out = 'obs_seq.processed', - first_obs_days = -1, - first_obs_seconds = -1, - last_obs_days = -1, - last_obs_seconds = -1, - print_only = .false., - gregorian_cal = .false. - / - -&obs_diag_nml - obs_sequence_name = 'obs_seq.final', - bin_width_days = -1, - bin_width_seconds = -1, - init_skip_days = 0, - init_skip_seconds = 0, - Nregions = 3, - trusted_obs = 'null', - lonlim1 = 0.00, 0.00, 0.50 - lonlim2 = 1.01, 0.50, 1.01 - reg_names = 'whole', 'yin', 'yang' - create_rank_histogram = .true., - outliers_in_histogram = .true., - use_zero_error_obs = .false., - verbose = .false. - / - -&state_vector_io_nml - / - -&model_mod_check_nml - input_state_files = 'input.nc' - output_state_files = 'mmc_output.nc' - test1thru = 0 - run_tests = 1,2,3,4,5,7 - x_ind = 42 - loc_of_interest = 0.3 - quantity_of_interest = 'QTY_STATE_VARIABLE' - interp_test_dx = 0.02 - interp_test_xrange = 0.0, 1.0 - verbose = .false. - / - -&quality_control_nml - input_qc_threshold = 3.0, - outlier_threshold = -1.0, -/ - -&location_nml - horiz_dist_only = .true. - approximate_distance = .false. - nlon = 71 - nlat = 36 - output_box_info = .true. -/ diff --git a/models/cice-scm2/work/quickbuild.sh b/models/cice-scm2/work/quickbuild.sh deleted file mode 100755 index e79b90dcb2..0000000000 --- a/models/cice-scm2/work/quickbuild.sh +++ /dev/null @@ -1,60 +0,0 @@ -#!/usr/bin/env bash - -# DART software - Copyright UCAR. This open source software is provided -# by UCAR, "as is", without charge, subject to all terms of use at -# http://www.image.ucar.edu/DAReS/DART/DART_download - -main() { - -export DART=$(git rev-parse --show-toplevel) -source "$DART"/build_templates/buildfunctions.sh - -MODEL=cice-scm2 -LOCATION=threed_sphere - - -programs=( -closest_member_tool -filter -model_mod_check -perfect_model_obs -) - -serial_programs=( -create_fixed_network_seq -create_obs_sequence -fill_inflation_restart -integrate_model -obs_common_subset -obs_diag -obs_sequence_tool -) - -model_programs=( -) - -model_serial_programs=( -dart_to_cice -) - -# quickbuild arguments -arguments "$@" - -# clean the directory -\rm -f -- *.o *.mod Makefile .cppdefs - -# build any NetCDF files from .cdl files -cdl_to_netcdf - -# build and run preprocess before making any other DART executables -buildpreprocess - -# build -buildit - -# clean up -\rm -f -- *.o *.mod - -} - -main "$@"