module mod_waves_only c --- Dany's/Jon's wave reading files are in here now [Tim: 2012.10.18] c --- So far there are c --- 1. wavewatch 3 from IFREMER (WAVES_WW3 flag) c --- 2. WAM North Sea from met.no (WAVES_WAMNSEA flag) c --- 3. old WAM (no flag) c --- swh, mwp, mwd are global variables that are now stored in here c --- also need to set waves_timeres (interval between wave readings) c --- WAVES_SKIP flag can also be used to skip readings #if defined (WAVES) #if defined (WAVES_ONLY) !#if defined (WAVES_TEST_READWAVES) ! use mod_common_ice ! !! - only need if want to check dfloe before it goes into ! !! waves_advect !#endif use mod_xc use mod_common_ice !use mod_common_wavesice!swh, mwp etc in there now implicit none ! real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: ! & ficem,hicem #if defined (WAVES_ICEOBS) integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv_conc,jpiv_conc real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist_conc real, parameter,private :: f_icemin = 0.05 !!check consistent with mod_common_wavesice.F !!should be a global var, but don't want to conflict with anything !!-should be the same as in m_icemodels_advect.F !! TODO put this into mod_common_ice.F !! but would need to be very careful #if defined (WAVES_IO_THICK1) !!SMOS pivot points integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv_thick,jpiv_thick real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist_thick #else/* Default = constant value */ !!use constant value !real, parameter,private :: h_fixed = 1.5 !!set ice thickness to constant value (chosen from small ice transect & drills, Fram Strait, Sept 2013) real, parameter,private :: h_fixed = 1.0 !!set ice thickness to constant value (test value) #endif #else/* ! defined WAVES_ICEOBS: use model daily averages */ type, private :: daily_averages_2 integer :: nrmem integer :: initday, inityear integer :: dumpday, dumpyear real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & fice,hice end type daily_averages_2 #endif character(len=*),parameter,private :: icedirec='./ICECONS_INPUT/' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_ice_thick_conc(dtime) implicit none real, intent(in) :: dtime #if defined (WAVES_ICEOBS) real, parameter :: fice_min = 0.15!! should this be the same as in !! m_icemodels_advect.F & mod_common_wavesice.F? !! TODO put this into mod_common_ice.F !! but would need to be very careful #if defined (WAVES_IO_THICK1) real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: htmp #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!read in conc #if defined (WAVES_IO_CONC2) call read_iceconc_osisaf(dtime) !get ficem from OSISAF #elif defined (WAVES_IO_CONC3) call read_iceconc_amsr2_3125(dtime) !get ficem from AMSR2 (3.125km) #else/* DEFAULT= AMSR-E (6.25km)*/ call read_iceconc_amsre6250(dtime) !get ficem from AMSR-E (6.25km) #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!read in thickness #if defined (WAVES_IO_THICK1) call read_icethick_smos(dtime) !get hicem from htmp = hicem ! SMOS thickness ! - NB this is effective thickness = conc*thickness hicem = 0.0 where (ficem.gt.fice_min) hicem = htmp/ficem #else/* DEFAULT = use fixed thickness */ hicem = 0.0 where (ficem.gt.fice_min) hicem = h_fixed #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #else/* ! defined WAVES_ICEOBS*/ call read_icecons_dailyave(dtime) #endif end subroutine read_ice_thick_conc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #if defined (WAVES_ICEOBS) #if defined (WAVES_IO_THICK1) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_pivots_smos !! lon-lat positions for concentrations from AMSR-E Bremen !use netcdf use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface !use mod_year_info, only: year_info, year_day,refyear implicit none !!pivot point files in SCRATCH/pivots !! - for nearest neighbour interpolation !! - these depend on the region !! - in preprocess.sh we link to the correct files character(len=*), parameter :: & fname_a = icedirec//'pivots/pivots_smos.a' & ,fname_b = icedirec//'pivots/pivots_smos.b' integer, parameter :: unit_piv = 926!!this can be anything I think character :: preambl(5)*79 character(len=*), parameter & :: data_name = 'SMOS thicknesses' integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv,jpiv real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist ! & ,ipivreal, jpivreal!!real values for output integer :: i,j logical :: pivf_exists include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!say if hycom is calling this routine: if (mnproc==1) then write(*,*) & "***************************************************" write(*,*) & "***** TW: Get pivot points for "//data_name//" ********" write(*,*) & "***************************************************" inquire(file=fname_a,exist=pivf_exists) print*,'TW: piv file = '//fname_a print*,'TW: piv file present?',pivf_exists endif ! use rdmonthch directly call zaiopf(fname_a, 'old',unit_piv) if (mnproc.eq.1) then ! .b file from 1st tile only open (unit=uoff+unit_piv,file=fname_b, & status='old', action='read') read (uoff+unit_piv,'(a79)') preambl endif !1st tile call preambl_print(preambl) call rdmonthck(util1, unit_piv,1) call rdmonthck(util2, unit_piv,2) call rdmonthck(util3, unit_piv,3) ! call xctilr( util1, 1,1, 1, 1, halo_ps) ! --- call rdmonth(util3, 990) if (mnproc.eq.1) then ! .b file from 1st tile only write(*,*)"ipiv file: ",fname_a close (unit=uoff+unit_piv) endif !1st tile call zaiocl(unit_piv) call xcsync(flush_lp) !! fill up the PP matrices !! - these are global variables and only need to be read once when !! initialising model (so doesn't depend on 'dtime') do j= 1,jj do i= 1,ii ipiv_thick(i,j) = util1(i,j) jpiv_thick(i,j) = util2(i,j) dist_thick(i,j) = util3(i,j) enddo !i enddo !j end subroutine read_pivots_smos !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_icethick_smos(dtime) !! lon-lat positions for concentrations from AMSR-E Bremen !use mod_xc !use mod_forcing_nersc !rdncrec use mod_hycom_nersc use netcdf use mod_xc ! HYCOM communication interface !use mod_za ! HYCOM I/O interface use mod_year_info, only: year_info, year_day,refyear implicit none real, intent(in) :: dtime character(len=*),parameter :: var_name = 'z' character(len=*),parameter & :: data_name = 'SMOS thicknesses' c c --- File where the lon's and lat's are located (in directory data_path): ! character(len=*),parameter :: lon_name = 'Longitudes' ! character(len=*),parameter :: lat_name = 'Latitudes' ! character(len=*), parameter :: ncname1 = ! & 'LongitudeLatitudeGrid-SMOS-Arctic.nc' character(len=80) :: ncfil1 c --- File where the data is kept: character(len=*), parameter :: data_path=icedirec//'SMOS/netcdf/' character(len=*), parameter :: & ncname2a = 'SMOS_Icethickness_north_', ncname2b = '.nc' character(len=8) :: cal_date!date eg 20110201 character(len=35) :: ncname2!eg SMOS_Icethickness_north_20110201.nc c --- Other temporary variables: integer :: nx2,ny2 integer :: ipiv,jpiv integer :: dimid, varid, ncid !!Declare longitude and latitude arrays !!-originals character(len=*), parameter :: ind_xname='x',ind_yname='y' integer, parameter :: ny = 896!no of y's integer, parameter :: nx = 608!no of x's real, dimension (nx,ny) :: & h0 !SMOS thicknesses ! & ,lon_h_obs !longitudes of the grid-points of the SMOS thicknesses ! & ,lat_h_obs !latitudes of the grid-points of the SMOS thicknesses integer :: i,j logical :: nan_test!,huge_test type(year_info) :: tt character(len=2) :: cdm_1is1!day in month (tt%cdm) starts from 0 ! - we want it to start from 1 include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!say if hycom is calling this routine: if (mnproc==1)write(*,*) & "***************************************************" if (mnproc==1)write(*,*) & "***** TW: Get lon-lat points for "//data_name//" ********" if (mnproc==1)write(*,*) & "***************************************************" !! thickness from observation !! set-up: get files and file/dim id's etc if (.false.) then!use example for testing cal_date = '20110201' else!get cal_date from dtime call year_day(dtime,refyear,tt,yrflag) write(cdm_1is1,'(i2.2)') tt%idm+1 cal_date = tt%cyy//tt%cmm//cdm_1is1 if (mnproc==1) print*,'TW: dtime ',dtime if (mnproc==1) print*,'TW: date ',cal_date end if ncname2 = ncname2a//cal_date//ncname2b ncfil1 = data_path//ncname2 if (mnproc==1)write(*,*)"SMOS path: ",data_path if (mnproc==1)write(*,*)"name of netcdf file (data): ",ncfil1 call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid)) call ncerr(nf90_inq_dimid(ncid,ind_xname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=nx2)) call ncerr(nf90_inq_dimid(ncid,ind_yname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=ny2)) if (nx/=nx2.or.ny/=ny2) then if (mnproc==1) then write(lp,'(a)') & 'TW: Dimension mismatch in data file (for '// & data_name//')' write(lp,'(a,2i5)') 'ny :',ny,ny2 write(lp,'(a,2i5)') 'nx :',nx,nx2 end if end if !! read in the variables; call ncerr(nf90_inq_varid(ncid,var_name,varid)) call ncerr(nf90_get_var (ncid,varid,h0)) call ncerr(nf90_close (ncid)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! hicem = 0.0 do j= 1,jj do i= 1,ii ipiv = ipiv_thick(i,j) jpiv = jpiv_thick(i,j) if (ipiv.gt.0) then !!in data grid hicem(i,j) = h0(ipiv,jpiv) end if !!missing values etc ('_FillValue'=NaN) !huge_test = (hicem(i,j).gt.100.0) nan_test = .not. & ( (hicem(i,j)+1.0).gt.(hicem(i,j)) ) if (nan_test) hicem(i,j) = 0.0 if ((i.eq.itest).and.(j.eq.jtest)) then print*, 'TW: checking SMOS reading' print*, 'TW: i,j = ',i,j print*, 'TW: ipiv,jpiv = ',ipiv,jpiv print*, 'TW: lon,lat = ',plon(i,j),plat(i,j) print*, 'TW: hicem = ',hicem(i,j) print*, 'TW: NaN? = ',nan_test !print*, 'TW: huge? = ',huge_test end if enddo !i enddo !j end subroutine read_icethick_smos !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif/*WAVES_IO_THICK1: SMOS thickness*/ #if defined (WAVES_IO_CONC2)/* OSISAF conc */ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_pivots_osisaf !! lon-lat positions for concentrations from AMSR-E Bremen !use mod_xc !use mod_forcing_nersc !rdncrec !use mod_hycom_nersc !use netcdf use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface !use mod_year_info, only: year_info, year_day,refyear implicit none !!pivot point files in SCRATCH/pivots !! - for nearest neighbour interpolation !! - these depend on the region !! - in preprocess.sh we link to the correct files character(len=*), parameter :: & fname_a = icedirec//'pivots/pivots_osisaf.a' & ,fname_b = icedirec//'pivots/pivots_osisaf.b' integer, parameter :: unit_piv = 926!!this can be anything I think character :: preambl(5)*79 character(len=*), parameter & :: data_name = 'OSISAF concentrations' integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv,jpiv real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist ! & ,ipivreal, jpivreal!!real values for output integer :: i,j logical :: pivf_exists include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!say if hycom is calling this routine: if (mnproc==1) then write(*,*) & "***************************************************" write(*,*) & "***** TW: Get pivot points for "//data_name//" ********" write(*,*) & "***************************************************" inquire(file=fname_a,exist=pivf_exists) print*,'TW: piv file = '//fname_a print*,'TW: piv file present?',pivf_exists endif ! use rdmonthch directly call zaiopf(fname_a, 'old',unit_piv) if (mnproc.eq.1) then ! .b file from 1st tile only open (unit=uoff+unit_piv,file=fname_b, & status='old', action='read') read (uoff+unit_piv,'(a79)') preambl endif !1st tile call preambl_print(preambl) call rdmonthck(util1, unit_piv,1) call rdmonthck(util2, unit_piv,2) call rdmonthck(util3, unit_piv,3) ! call xctilr( util1, 1,1, 1, 1, halo_ps) ! --- call rdmonth(util3, 990) if (mnproc.eq.1) then ! .b file from 1st tile only write(*,*)"ipiv file: ",fname_a close (unit=uoff+unit_piv) endif !1st tile call zaiocl(unit_piv) call xcsync(flush_lp) if (mnproc==1) print*,'TW: OSISAF pivot file read successfully' !! fill up the PP matrices !! - these are global variables and only need to be read once when !! initialising model (so doesn't depend on 'dtime') do j= 1,jj do i= 1,ii ipiv_conc(i,j) = util1(i,j) jpiv_conc(i,j) = util2(i,j) dist_conc(i,j) = util3(i,j) enddo !i enddo !j end subroutine read_pivots_osisaf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_iceconc_osisaf(dtime) !! lon-lat positions for concentrations from AMSR-E Bremen !use mod_xc !use mod_forcing_nersc !rdncrec use mod_hycom_nersc use netcdf use mod_xc ! HYCOM communication interface !use mod_za ! HYCOM I/O interface use mod_year_info, only: year_info, year_day,refyear implicit none real, intent(in) :: dtime character(len=*),parameter :: var_name = 'ice_conc' !!data has a confidence level associated with it character(len=*),parameter :: con_name = 'confidence_level' real,parameter :: con_lim = 3.0 character(len=*), parameter & :: data_name = 'OSISAF concentrations' c c --- File where the data is kept: character(len=*), parameter :: & data_path=icedirec//'/OSISAF/' character(len=80) :: ncfil1 character(len=*), parameter :: & ncname2a = 'ice_conc_nh_polstere-100_multi_' character(len=*), parameter :: ncname2b = '1200.nc' character(len=8) :: cal_date!date eg 20130901 character(len=46) :: ncname2!eg ice_conc_nh_polstere-100_multi_201309011200.nc c --- Other temporary variables: integer :: nx2,ny2 integer :: dimid, varid, ncid integer :: ipiv,jpiv !!Declare longitude and latitude arrays !!-originals character(len=*), parameter :: ind_xname='xc',ind_yname='yc' integer, parameter :: ny = 1120!no of y's integer, parameter :: nx = 760!no of x's real, dimension (nx,ny) :: & conc_obs0,var_con !OSISAF concentrations and confidence levels integer :: i,j logical :: badval type(year_info) :: tt character(len=2) :: cdm_1is1!day in month (tt%cdm) starts from 0 ! - we want it to start from 1 include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! concentration from observation !! set-up: get files and file/dim id's etc if (.false.) then!use example for testing cal_date = '20110201' else!get cal_date from dtime call year_day(dtime,refyear,tt,yrflag) write(cdm_1is1,'(i2.2)') tt%idm+1 !!1st of month tt%cdm is 0, !! but tt%idm is 1 cal_date = tt%cyy//tt%cmm//cdm_1is1 if (mnproc==1) print*,'TW: dtime ',dtime if (mnproc==1) print*,'TW: date ',cal_date end if ncname2 = ncname2a//cal_date//ncname2b ncfil1 = data_path//ncname2 if (mnproc==1)write(*,*)"name of netcdf file (data): ",ncfil1 call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))!!**OPENED call ncerr(nf90_inq_dimid(ncid,ind_xname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=nx2)) call ncerr(nf90_inq_dimid(ncid,ind_yname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=ny2)) if (nx/=nx2.or.ny/=ny2) then if (mnproc==1) then write(lp,'(a)') & 'TW: Dimension mismatch in data file (for '// & data_name//')' write(lp,'(a,2i5)') 'ny :',ny,ny2 write(lp,'(a,2i5)') 'nx :',nx,nx2 end if end if !! read in the conc; call ncerr(nf90_inq_varid(ncid,var_name,varid)) call ncerr(nf90_get_var (ncid,varid,conc_obs0)) !! read in the conf level; call ncerr(nf90_inq_varid(ncid,con_name,varid)) call ncerr(nf90_get_var (ncid,varid,var_con)) call ncerr(nf90_close (ncid))!!**CLOSED call rdncrec(ncfil1,var_name,conc_obs0,nx,ny,1)!applies scale factor etc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ficem = 0.0 do j= 1,jj do i= 1,ii ipiv = ipiv_conc(i,j) jpiv = jpiv_conc(i,j) if (ipiv.gt.0) then ficem(i,j) = conc_obs0(ipiv,jpiv) end if !!missing values etc badval = ((ficem(i,j).lt.0.0).or. & (var_con(ipiv,jpiv).lt.3)) if ((i==itest).and.(j==jtest)) then print*,'TWijt: i,j:',i+i0,j+j0 print*,'TWijt: ipiv,jpiv:',i+i0,j+j0 print*,'TWijt: conc,conc0,conf:',ficem(i,j), & conc_obs0(ipiv,jpiv),var_con(ipiv,jpiv) print*,'TWijt: badval:',badval end if if (badval) ficem(i,j) = 0.0 enddo !i enddo !j !!scale to 0<=ficem<=1 ficem = ficem/100.0 end subroutine read_iceconc_osisaf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #elif defined (WAVES_IO_CONC3)/* AMSR2 (3.125km) conc */ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_pivots_amsr2_3125 !! lon-lat positions for concentrations from AMSR2 (3.125km) !use mod_xc !use mod_forcing_nersc !rdncrec !use mod_hycom_nersc !use netcdf use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface !use mod_year_info, only: year_info, year_day,refyear implicit none !!pivot point files in SCRATCH/pivots !! - for nearest neighbour interpolation !! - these depend on the region !! - in preprocess.sh we link to the correct files character(len=*), parameter :: & fname_a = icedirec//'pivots/pivots_amsr2_3125.a' & ,fname_b = icedirec//'pivots/pivots_amsr2_3125.b' integer, parameter :: unit_piv = 926!!this can be anything I think character :: preambl(5)*79 character(len=*), parameter & :: data_name = 'AMSR2 (3.125km) concentrations' integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv,jpiv real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist ! & ,ipivreal, jpivreal!!real values for output integer :: i,j logical :: pivf_exists include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!say if hycom is calling this routine: if (mnproc==1) then write(*,*) & "***************************************************" write(*,*) & "***** TW: Get pivot points for "//data_name//" ********" write(*,*) & "***************************************************" inquire(file=fname_a,exist=pivf_exists) print*,'TW: piv file = '//fname_a print*,'TW: piv file present?',pivf_exists endif ! use rdmonthch directly call zaiopf(fname_a, 'old',unit_piv) if (mnproc.eq.1) then ! .b file from 1st tile only open (unit=uoff+unit_piv,file=fname_b, & status='old', action='read') read (uoff+unit_piv,'(a79)') preambl endif !1st tile call preambl_print(preambl) call rdmonthck(util1, unit_piv,1) call rdmonthck(util2, unit_piv,2) call rdmonthck(util3, unit_piv,3) ! call xctilr( util1, 1,1, 1, 1, halo_ps) ! --- call rdmonth(util3, 990) if (mnproc.eq.1) then ! .b file from 1st tile only write(*,*)"ipiv file: ",fname_a close (unit=uoff+unit_piv) endif !1st tile call zaiocl(unit_piv) call xcsync(flush_lp) if (mnproc==1) print*,'TW: ',data_name, & ' pivot file read successfully' !! fill up the PP matrices !! - these are global variables and only need to be read once when !! initialising model (so doesn't depend on 'dtime') do j= 1,jj do i= 1,ii ipiv_conc(i,j) = util1(i,j) jpiv_conc(i,j) = util2(i,j) dist_conc(i,j) = util3(i,j) enddo !i enddo !j end subroutine read_pivots_amsr2_3125 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_iceconc_amsr2_3125(dtime) !! Concentrations from AMSR2 Bremen (3.125km) !use mod_xc !use mod_forcing_nersc !rdncrec use mod_hycom_nersc use netcdf use mod_xc ! HYCOM communication interface !use mod_za ! HYCOM I/O interface use mod_year_info, only: year_info, year_day,refyear implicit none real, intent(in) :: dtime character(len=*),parameter :: var_name = 'sea_ice_concentration' ! !!data has a confidence level associated with it ! character(len=*),parameter :: con_name = 'confidence_level' ! real,parameter :: con_lim = 3.0 character(len=*), parameter & :: data_name = 'AMSR2 (3.125km) concentrations' c c --- File where the data is kept: character(len=*), parameter :: & data_path=icedirec//'/AMSR2_3125/' character(len=80) :: ncfil1 character(len=*), parameter :: & ncname2a = 'Arc_' character(len=*), parameter :: ncname2b = '_res3.125_pyres.nc' character(len=8) :: cal_date!date eg 20130901 character(len=46) :: ncname2!eg Arc_20130921_res3.125_pyres.nc c --- Other temporary variables: integer :: nx2,ny2 integer :: dimid, varid, ncid integer :: ipiv,jpiv !!Declare longitude and latitude arrays !!-originals character(len=*), parameter :: ind_xname='x',ind_yname='y' integer, parameter :: ny = 3584!no of y's integer, parameter :: nx = 2432!no of x's real, dimension (nx,ny) :: & conc_obs0 !AMSR2 (3.125km) concentrations integer :: i,j logical :: badval type(year_info) :: tt character(len=2) :: cdm_1is1!day in month (tt%cdm) starts from 0 ! - we want it to start from 1 include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! concentration from observation !! set-up: get files and file/dim id's etc if (.false.) then!use example for testing cal_date = '20110201' else!get cal_date from dtime call year_day(dtime,refyear,tt,yrflag) write(cdm_1is1,'(i2.2)') tt%idm+1 !!1st of month tt%cdm is 0, !! but tt%idm is 1 cal_date = tt%cyy//tt%cmm//cdm_1is1 if (mnproc==1) print*,'TW: dtime ',dtime if (mnproc==1) print*,'TW: date ',cal_date end if ncname2 = ncname2a//cal_date//ncname2b ncfil1 = data_path//'Arc_'//tt%cyy//'/'//ncname2 if (mnproc==1)write(*,*)"name of netcdf file (data): ",ncfil1 call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))!!**OPENED call ncerr(nf90_inq_dimid(ncid,ind_xname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=nx2)) call ncerr(nf90_inq_dimid(ncid,ind_yname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=ny2)) if (nx/=nx2.or.ny/=ny2) then if (mnproc==1) then write(lp,'(a)') & 'TW: Dimension mismatch in data file (for '// & data_name//')' write(lp,'(a,2i5)') 'ny :',ny,ny2 write(lp,'(a,2i5)') 'nx :',nx,nx2 end if end if !! read in the conc; call ncerr(nf90_inq_varid(ncid,var_name,varid)) call ncerr(nf90_get_var (ncid,varid,conc_obs0)) ! !! read in the conf level; ! call ncerr(nf90_inq_varid(ncid,con_name,varid)) ! call ncerr(nf90_get_var (ncid,varid,var_con)) ! call ncerr(nf90_close (ncid))!!**CLOSED call rdncrec(ncfil1,var_name,conc_obs0,nx,ny,1)!applies scale factor etc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ficem = 0.0 do j= 1,jj do i= 1,ii ipiv = ipiv_conc(i,j) jpiv = jpiv_conc(i,j) if (ipiv.gt.0) then ficem(i,j) = conc_obs0(ipiv,jpiv) end if !!missing values etc badval = ((ficem(i,j).lt.0.0).or. & (ficem(i,j).gt.100.0)) if ((i==itest).and.(j==jtest)) then print*,'TWijt: i,j:',i+i0,j+j0 print*,'TWijt: ipiv,jpiv:',i+i0,j+j0 print*,'TWijt: conc,conc0,conf:',ficem(i,j), & conc_obs0(ipiv,jpiv)!,var_con(ipiv,jpiv) print*,'TWijt: badval:',badval end if if (badval) ficem(i,j) = 0.0 enddo !i enddo !j !!scale to 0<=ficem<=1 ficem = ficem/100.0 end subroutine read_iceconc_amsr2_3125 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #else/* DEFAULT = AMSR-E (6.25km)*/ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_pivots_amsre6250 !! lon-lat positions for concentrations from AMSR-E Bremen !use mod_xc !use mod_forcing_nersc !rdncrec !use mod_hycom_nersc !use netcdf use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface !use mod_year_info, only: year_info, year_day,refyear implicit none !!pivot point files in SCRATCH/pivots !! - for nearest neighbour interpolation !! - these depend on the region !! - in preprocess.sh we link to the correct files character(len=*), parameter :: & fname_a = icedirec//'pivots/pivots_amsre.a' & ,fname_b = icedirec//'pivots/pivots_amsre.b' integer, parameter :: unit_piv = 926!!this can be anything I think character :: preambl(5)*79 character(len=*), parameter & :: data_name = 'AMSR-E concentrations' integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ipiv,jpiv real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dist ! & ,ipivreal, jpivreal!!real values for output integer :: i,j logical :: pivf_exists include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!say if hycom is calling this routine: if (mnproc==1) then write(*,*) & "***************************************************" write(*,*) & "***** TW: Get pivot points for "//data_name//" ********" write(*,*) & "***************************************************" inquire(file=fname_a,exist=pivf_exists) print*,'TW: piv file = '//fname_a print*,'TW: piv file present?',pivf_exists endif ! use rdmonthch directly call zaiopf(fname_a, 'old',unit_piv) if (mnproc.eq.1) then ! .b file from 1st tile only open (unit=uoff+unit_piv,file=fname_b, & status='old', action='read') read (uoff+unit_piv,'(a79)') preambl endif !1st tile call preambl_print(preambl) call rdmonthck(util1, unit_piv,1) call rdmonthck(util2, unit_piv,2) call rdmonthck(util3, unit_piv,3) ! call xctilr( util1, 1,1, 1, 1, halo_ps) ! --- call rdmonth(util3, 990) if (mnproc.eq.1) then ! .b file from 1st tile only write(*,*)"ipiv file: ",fname_a close (unit=uoff+unit_piv) endif !1st tile call zaiocl(unit_piv) call xcsync(flush_lp) if (mnproc==1) print*,'TW: AMSR-E pivot file read successfully' !! fill up the PP matrices !! - these are global variables and only need to be read once when !! initialising model (so doesn't depend on 'dtime') do j= 1,jj do i= 1,ii ipiv_conc(i,j) = util1(i,j) jpiv_conc(i,j) = util2(i,j) dist_conc(i,j) = util3(i,j) enddo !i enddo !j end subroutine read_pivots_amsre6250 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_iceconc_amsre6250(dtime) !! lon-lat positions for concentrations from AMSR-E Bremen !use mod_xc !use mod_forcing_nersc !rdncrec use mod_hycom_nersc use netcdf use mod_xc ! HYCOM communication interface !use mod_za ! HYCOM I/O interface use mod_year_info, only: year_info, year_day,refyear implicit none real, intent(in) :: dtime character(len=*),parameter :: var_name = 'sea_ice_concentration' character(len=*), parameter & :: data_name = 'AMSR-E concentrations' c c --- File where the lon's and lat's are located (in directory data_path): !! !is this needed? !! character(len=*), parameter :: ncname1 = !! & 'LongitudeLatitudeGrid-n6250-Arctic.nc' !! character(len=*),parameter :: lon_name = 'Longitudes' !! character(len=*),parameter :: lat_name = 'Latitudes' c --- File where the data is kept: character(len=*), parameter :: & data_path=icedirec//'/AMSRE6.25km/netcdf/' character(len=80) :: ncfil1 character(len=*), parameter :: ncname2a = 'asi-n6250-' character(len=*), parameter :: ncname2b = '-v5i.nc' character(len=8) :: cal_date!date eg 20110201 character(len=25) :: ncname2!eg asi-n6250-20110201-v5i.nc c --- Other temporary variables: integer :: nx2,ny2 integer :: dimid, varid, ncid integer :: ipiv,jpiv !!Declare longitude and latitude arrays !!-originals character(len=*), parameter :: ind_xname='x',ind_yname='y' integer, parameter :: ny = 1792!no of y's integer, parameter :: nx = 1216!no of x's real, dimension (nx,ny) :: & conc_obs0 !AMSR-E concentrations !! & ,lon_conc_obs !longitudes of the grid-points of the AMSR-E concentrations !! & ,lat_conc_obs !latitudes of the grid-points of the AMSR-E concentrations integer :: i,j type(year_info) :: tt character(len=2) :: cdm_1is1!day in month (tt%cdm) starts from 0 ! - we want it to start from 1 include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! concentration from observation !! set-up: get files and file/dim id's etc if (.false.) then!use example for testing cal_date = '20110201' else!get cal_date from dtime call year_day(dtime,refyear,tt,yrflag) write(cdm_1is1,'(i2.2)') tt%idm+1 !!1st of month tt%cdm is 0, !! but tt%idm is 1 cal_date = tt%cyy//tt%cmm//cdm_1is1 if (mnproc==1) print*,'TW: dtime ',dtime if (mnproc==1) print*,'TW: date ',cal_date end if ncname2 = ncname2a//cal_date//ncname2b ncfil1 = data_path//ncname2 if (mnproc==1)write(*,*)"name of netcdf file (data): ",ncfil1 call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid)) call ncerr(nf90_inq_dimid(ncid,ind_xname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=nx2)) call ncerr(nf90_inq_dimid(ncid,ind_yname,dimid)) call ncerr(nf90_inquire_dimension(ncid,dimid,len=ny2)) if (nx/=nx2.or.ny/=ny2) then if (mnproc==1) then write(lp,'(a)') & 'TW: Dimension mismatch in data file (for '// & data_name//')' write(lp,'(a,2i5)') 'ny :',ny,ny2 write(lp,'(a,2i5)') 'nx :',nx,nx2 end if end if !! read in the variables; call ncerr(nf90_inq_varid(ncid,var_name,varid)) call ncerr(nf90_get_var (ncid,varid,conc_obs0)) call ncerr(nf90_close (ncid)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ficem = 0.0 do j= 1,jj do i= 1,ii ipiv = ipiv_conc(i,j) jpiv = jpiv_conc(i,j) if (ipiv.gt.0) then ficem(i,j) = conc_obs0(ipiv,jpiv) end if !!missing values etc if (ficem(i,j).gt.100.0) ficem(i,j) = 0.0 enddo !i enddo !j !!scale to 0<=ficem<=1 ficem = ficem/100.0 end subroutine read_iceconc_amsre6250 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif/* CONC DATA OPTIONS */ #else/* ! defined WAVES_ICEOBS: use model daily averages */ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_icecons_dailyave(dtime) use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface !use mod_daily_average use mod_year_info, only: year_info, year_day,refyear implicit none type(daily_averages_2) :: da real, intent(in) :: dtime character(len=80) :: fname_a,fname_b character(len=8) :: cal_date character(len=*),parameter :: & fname_start = 'dailyave_waves_input_' & ,data_path = './ICECONS_INPUT/DAILYAVE_WAVES_INPUT/' integer, parameter :: unit_piv = 926,irec_f=10,irec_h=11!!unit_piv can be anything I think character :: preambl(1016)*80 character(len=*), parameter :: & data_name = 'daily average conc & thickness' ! integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: ! & ipiv,jpiv ! real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: ! & dist ! & ,ipivreal, jpivreal!!real values for output integer :: i,j logical :: pivf_exists type(year_info) :: tt character(len=2) :: cdm_1is1!day in month (tt%cdm) starts from 0 ! - we want it to start from 1 include 'common_blocks.h' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (.false.) then cal_date = '20100624' else!work out cal_date from dtime call year_day(dtime,refyear,tt,yrflag) write(cdm_1is1,'(i2.2)') tt%idm+1 !!1st of month tt%cdm is 0, !! but tt%idm is 1 cal_date = tt%cyy//tt%cmm//cdm_1is1 if (mnproc==1) print*,'TW: dtime ',dtime if (mnproc==1) print*,'TW: date ',cal_date end if fname_a = data_path//fname_start//cal_date//'.a' fname_b = data_path//fname_start//cal_date//'.b' !!say if hycom is calling this routine: if (mnproc==1) then write(*,*) & "***************************************************" write(*,*) & "***** TW: Read "//data_name//" ********" inquire(file=fname_a,exist=pivf_exists) print*,'TW: Daily ave input file = '//fname_a print*,'TW: Daily ave input file present?',pivf_exists write(*,*) & "***************************************************" endif call daily_ave_read_2(da,fname_a,fname_b) ficem = da%fice hicem = da%hice end subroutine read_icecons_dailyave !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine daily_ave_read_2(datmp, filenamea,filenameb, & rmem,rstep,rtime) use mod_za implicit none type(daily_averages_2), intent(inout) :: datmp !integer, intent(in) :: inityear, initday,dumpyear,dumpday !TW integer, optional, intent(in) :: rmem integer, optional, intent(out) :: rstep real, optional, intent(out) :: rtime real :: xmin,xmax,xmin2,xmax2,coord,loctime integer :: idummy,jdummy,l_itdm,l_jtdm,l_kdm,k,nop, ios,n, & locnstep, lociversn, lociexpt, locyrflag, & ktr,loop_no character(len=80) :: locctitle(4) !character(len=60) :: filenamea,filenameb, filebase !TW character(len=80), intent(in) :: filenamea,filenameb character(len=8) :: char8 logical :: !no2dfld,no3dfld & got_fice,got_hice,got_ice_fh,no_read real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & dummy_fld include 'common_blocks.h' c c --- Set filenames !TW ! filebase=dailyfilebase(inityear, initday, ! & dumpyear, dumpday, rmem) ! filenamea=trim(filebase)//'.a' ! filenameb=trim(filebase)//'.b' if (mnproc==1) print *, & 'TW: a,b files='//trim(filenamea)//','//trim(filenameb) !TW if (mnproc==1) then write(lp,'(a)') 'Reading '// & trim(filenamea)//', '//trim(filenameb) !TW ! & trim(filebase)//'.[ab]' !TW call flush(lp) end if c c --- Open files nop=710 open (unit=nop,file=trim(filenameb), & status='old') call zaiopf(trim(filenamea),'old',nop) c c --- Read .b Header read(nop,116) locctitle,lociversn,lociexpt,locyrflag, & l_itdm,l_jtdm, & l_kdm,datmp%inityear,datmp%initday, & datmp%dumpyear,datmp%dumpday,datmp%nrmem c c -- Security Check if (l_itdm/=itdm.or.l_jtdm/=jtdm.or.l_kdm/=kdm) then Print *,'Grid dimensions are wrong!!!' call xcstop ('(mod_daily_average:daily_ave_read)') stop '(mod_daily_average:daily_ave_read)' end if c c --- Read the local ensemble members into ".b" header do n=1,1000 read(nop,118) idummy,jdummy!datmp%lmem(n) end do read(nop,119) if (mnproc==1) print*,'TW (daily_ave_read_2): made it to here (2)' c c --- Start reading fields Read until we fail ios=0 loop_no = 0 got_fice = .false. got_hice = .false. got_ice_fh = .false. do while ((ios==0).and.(.not.got_ice_fh)) loop_no = loop_no+1 no_read = .true. if (mnproc==1) print*, 'TW: ios loops',loop_no read (nop,117,iostat=ios) char8,locnstep,loctime, & k,coord,xmin2,xmax2 if (mnproc==1) print*, 'TW: char8,k,ios = ',char8,k,ios if (present(rstep)) then rstep = locnstep end if if (present(rtime)) then rtime = loctime end if ! no2dfld=.false. ! no3dfld=.false. ktr=0 if (mnproc==1) print*, & 'TW (daily_ave_read_2): made it to here (2b)' if (mnproc==1) print*,'TW: ios, k: ',ios, k if (ios==0) then c --- Read 2D fields if (k==0) then if (adjustl(trim(char8))=='ubavg') then if (mnproc==1) print*, & 'TW (daily_ave_read_2): made it to here (2c)' call zaiord(dummy_fld(1-nbdy,1-nbdy),iu,.false., & xmin,xmax, nop) if (mnproc==1) print*, & 'TW (daily_ave_read_2): made it to here (2d)' no_read = .false. else if (adjustl(trim(char8))=='vbavg') then call zaiord(dummy_fld(1-nbdy,1-nbdy),iv,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='ssh' ) then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. !#if defined (DAILY_ENSVAR) ! else if (adjustl(trim(char8))=='ssh_sq' ) then ! call zaiord(datmp%ssh_sq(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='fice_sq' ) then ! call zaiord(datmp%fice_sq(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='hice_sq' ) then ! call zaiord(datmp%hice_sq(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#endif else if (adjustl(trim(char8))=='surflx') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='swflx') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='salflx') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='emnp') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='taux') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. else if (adjustl(trim(char8))=='tauy') then call zaiord(dummy_fld(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) no_read = .false. !#if defined (ICE) || defined(ICESTATE) else if (adjustl(trim(char8))=='fice') then ! if (adjustl(trim(char8))=='fice') then call zaiord(datmp%fice(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) got_fice = .true. no_read = .false. else if (adjustl(trim(char8))=='hice') then call zaiord(datmp%hice(1-nbdy,1-nbdy),ip,.false., & xmin,xmax, nop) got_hice = .true. no_read = .false. ! else if (adjustl(trim(char8))=='hsnw') then ! call zaiord(datmp%hsnw(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='surflxs') then ! call zaiord(datmp%srfheatflux(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='albedo') then ! call zaiord(datmp%srfalb(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='htndncy') then ! call zaiord(datmp%htndncy(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#if defined(ALBSNW_EVOL) ! else if (adjustl(trim(char8))=='fsnw_cov') then ! call zaiord(datmp%srfsnw(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#endif !#endif !#if defined(EVP) ! else if (adjustl(trim(char8))=='uice') then ! call zaiord(datmp%uice(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='vice') then ! call zaiord(datmp%vice(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='spdice') then ! call zaiord(datmp%wice(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='tauxi ') then ! call zaiord(datmp%tauxice(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='tauyi ') then ! call zaiord(datmp%tauyice(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#if defined(ICE_DYN_DIAG) ! else if (adjustl(trim(char8))=='sigp ') then ! call zaiord(datmp%sigp(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='sigm ') then ! call zaiord(datmp%sigm(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='sig12 ') then ! call zaiord(datmp%sig12(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='pr ') then ! call zaiord(datmp%pr(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='strI ') then ! call zaiord(datmp%strI(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='strII ') then ! call zaiord(datmp%strII(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#endif !#if defined(TEST_ICE_AGE) ! else if (adjustl(trim(char8))=='fy_frac') then ! call zaiord(datmp%fy_frac(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='fy_age ') then ! call zaiord(datmp%fy_age(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='rdg_frac') then ! call zaiord(datmp%rdg_frac(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! !#endif !#if defined(PARAM_EST) ! else if (adjustl(trim(char8))=='sstb') then ! call zaiord(datmp%sstb(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='msshb ') then ! call zaiord(datmp%msshb(1-nbdy,1-nbdy),ip,.false., ! & xmin,xmax, nop) !#endif !#endif ! else ! no2dfld=.true. end if c --- Read 3D fields ! else if (k<=kdm .and. k>0) then ! ! if (adjustl(trim(char8))=='utot') then ! call zaiord(datmp%u(1-nbdy,1-nbdy,k),iu,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='vtot') then ! call zaiord(datmp%v(1-nbdy,1-nbdy,k),iv,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='saln' ) then ! call zaiord(datmp%s(1-nbdy,1-nbdy,k),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='temp') then ! call zaiord(datmp%t(1-nbdy,1-nbdy,k),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='pres') then ! call zaiord(datmp%p(1-nbdy,1-nbdy,k),ip,.false., ! & xmin,xmax, nop) ! else if (adjustl(trim(char8))=='kinetic') then ! call zaiord(datmp%w(1-nbdy,1-nbdy,k),ip,.false., ! & xmin,xmax, nop) !#if defined (NOR05) ! elseif (adjustl(trim(char8))=='NIT') then ! ktr=init ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='PHO') then ! ktr=ipho ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='SIL') then ! ktr=isil ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='DET') then ! ktr=idet ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='SIS') then ! ktr=isis ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='FLA') then ! ktr=ifla ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='DIA') then ! ktr=idia ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='OXY') then ! ktr=ioxy ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='SED') then ! ktr=ised ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='YEL') then ! ktr=iyel ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='DETP') then ! ktr=idetp ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='MICRO') then ! ktr=imicro ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='MESO') then ! ktr=imeso ! call zaiord(datmp%trcr_ave(1-nbdy,1-nbdy,k,ktr),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='primprod') then ! call zaiord(datmp%grosspp_ave(1-nbdy,1-nbdy,k),ip, ! & .false.,xmin,xmax, nop) ! elseif (adjustl(trim(char8))=='netpp') then ! call zaiord(datmp%netpp_ave(1-nbdy,1-nbdy,k),ip, ! & .false.,xmin,xmax, nop) !#endif /* NOR05 */ ! else ! no3dfld=.true. ! end if !c --- We just read a field with k> kdm or k<0. No good ! else ! if (mnproc==1) then ! write(lp,'(a)') 'Illegal specification of k' ! write(lp,'(a,i3)') 'Offending field '//char8,k ! end if ! call xcstop('(mod_daily_average:daily_ave_read)') ! stop '(mod_daily_average:daily_ave_read)' end if!k==0 c c --- Sanity check: ! if (no2dfld.and.no3dfld) then ! if (mnproc==1) then ! write(lp,'(a)') 'No flds corresponding to ', ! & adjustl(trim(char8)), ' found' ! end if ! call xcstop('(mod_daily_average:daily_ave_read)') ! stop '(mod_daily_average:daily_ave_read)' ! end if c if (.not.no_read) then c --- Check for coherency between min and maxes from a and b files if (abs(xmin-xmin2)>abs(xmin)*1e-4 .or. & abs(xmax-xmax2)>abs(xmax)*1e-4) then if (mnproc.eq.1) then write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') & 'error - .a and .b files not consistent:', & '.a,.b min = ',xmin,xmin2 ,xmin2-xmin , & '.a,.b max = ',xmax,xmax2 ,xmax2-xmax print *,'Variable :',char8,k endif call xcstop('(mod_daily_average:daily_ave_read)') stop '(mod_daily_average:daily_ave_read)' end if end if end if!ios==0 got_ice_fh = (got_fice.and.got_hice) if (mnproc==1) print*,'TW (DAR_2):', & got_fice,got_hice,got_ice_fh enddo!while ios==0 if (mnproc==1) print*,'TW (daily_ave_read_2): made it to here (3)' c c --- Finished - close files call zaiocl(nop) if (mnproc==1) close(nop) if (mnproc==1) print*,'TW (daily_ave_read_2): made it to here (4)' c 116 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & i5,4x,'''kdm '' = Vertical array size'/ & i5,4x,'''syear '' = Year of integration start '/ & i5,4x,'''sday '' = Day of integration start'/ & i5,4x,'''dyear '' = Year of this dump '/ & i5,4x,'''dday '' = Day of this dump '/ & i5,4x,'''count '' = Ensemble counter ') C & 'field time step model day', C & ' k dens min max') 117 format (a8,' = ',i11,f11.2,i3,f7.3,1p2e16.7) 118 format ('member ',i5.5,' = ',i1,' Ensemble member flag') 119 format('field time step model day', & ' k dens min max') if (mnproc==1) print*,'TW (daily_ave_read_2): made it to here (5)' end subroutine daily_ave_read_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif/* ! WAVES_ICEOBS */ #endif/* WAVES_ONLY */ #endif/* WAVES */ end module mod_waves_only