! source file: /sfs/fs1/work-geomar6/smomw113/UVic_ESCM/2.9_dk_0814/updates/07_aouassim/source/common/sstdata_panik.F subroutine sstdata_panik(locali,localj,sstreturn) !======================================================================= ! routine to read and interpolate three dimensional SST data for ASE experiments ! wkoeve@geomar.de, 22.06.16, experimental code so far, under development ! This routine shall be called from gasbc.F, similar to cfcdata.f. !======================================================================= implicit none character(120) :: fname, name, new_file_name, text integer iou, n, ln, ib(10), ic(10) integer locali, localj, c,i,j integer :: AllocateStatus, DeAllocateStatus logical inqvardef, exists real data_time, tim(3), wt1, wt3 real sstreturn real, allocatable :: time(:),tmpijkl(:,:,:,:),datasst(:,:,:) ! All variable values stored via save are available each time sstdata() ! is being entered, kind of internal global variables. include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "calendar.h" include "cembm.h" include "switch.h" include "tmngr.h" real dat(imtm2,jmtm2,3) save dat, datasst, ln, tim, time if (.not. allocated (time)) then ! variable time has not yet been alocated, hence this section of the code is only entered ! once during the model run name = "O_sstase.nc" fname = new_file_name (name) inquire (file=trim(fname), exist=exists) if (.not. exists) then print*, "==> Fatal Error: ", trim(fname), " does not exist." print*, "stopping UVic in sstdata.F" ! There shall be no fallback optipon: if the input file does not exist ! we simply stop UVic. stop else call openfile (fname, iou) call getdimlen ('MYTIME', iou, ln) write(20,*) 'sstdata: dim of Mytime', ln allocate ( time(ln) ) ! allocate ( tmpijkl(ln,1,jmtm2,imtm2), STAT = AllocateStatus ) allocate ( tmpijkl(imtm2,jmtm2,1,ln), STAT = AllocateStatus ) IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" allocate ( datasst(imtm2,jmtm2,ln), STAT = AllocateStatus ) IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" ! tmpijkl(1:ln,1,1:jmtm2,1:imtm2) = c0 tmpijkl(1:imtm2,1:jmtm2,1,1:ln) = c0 datasst(1:imtm2,1:jmtm2,1:ln) = c0 ib(:) = 1 ic(:) = ln ! ic(1) = ln ! ic(2) = 1 ! ic(3) = jmtm2 ! ic(4) = imtm2 ic(4) = ln ic(3) = 1 ic(2) = jmtm2 ic(1) = imtm2 ! note that so far the file has been written with ferret, ! i.e. with an offset of 1700 yrs call getvara ('MYTIME', iou, ln, 1, ln, time, c1, c0) text = 'years' call getatttext (iou, 'MYTIME', 'units', text) if (trim(text) .eq. "days since 1-1-1") & time(:) = time(:)/yrlen - 1. if (trim(text) .eq. "days since 0-1-1") & time(:) = time(:)/yrlen if (trim(text) .eq. "years since 1-1-1") & time(:) = time(:) - 1. do c=1,ln time(c)=time(c)+1700 ! write(20,*) time(c) enddo exists = inqvardef('O_SSTFORC', iou) if (.not. exists) then print*, "==> Fatal Error: O_sstforc data does not exist." stop else call getvara ('O_SSTFORC', iou, imtm2*jmtm2*ln, ib, ic &, tmpijkl, c1, c0) do i=1,imtm2 do j=1,jmtm2 do c=1,ln datasst(i,j,c) = ! & tmpijkl(c,1,j,i) & tmpijkl(i,j,1,c) enddo enddo enddo ! DEALLOCATE (tmpijkl, STAT = DeAllocateStatus) write(20,*) & 'Note: O_sstforc data have been loaded into memory.' endif endif tim(:) = time(1) do j=1,jmtm2 write(20,*) 'j = ',j, ' -------------------------' do i=1,imtm2 write(20,'(F7.1)',advance="no") datasst(i,j,1) dat(i,j,:) = datasst(i,j,1) enddo enddo endif data_time = year0 + accel_yr0 + (relyr - accel_yr0)*accel tim(2) = min(time(ln), max(time(1), data_time)) ! The following part is the time-interpolation part. ! For time values outside the range of the sstdata-file, I leave it ! as in the cfc case. Note, that this causes the sst-for-ase to be ! not very well defined outside this range with respect to seasonal ! SST changes! if (tim(2) .le. time(1)) then dat(locali,localj,2) = datasst(locali,localj,1) elseif (tim(2) .ge. time(ln)) then dat(locali,localj,2) = datasst(locali,localj,ln) else if (tim(2) .gt. tim(3)) then do n=2,ln if (time(n-1) .le. tim(2) .and. time(n) .ge. tim(2)) then tim(1) = time(n-1) dat(locali,localj,1) = datasst(locali,localj,n-1) tim(3) = time(n) dat(locali,localj,3) = datasst(locali,localj,n) endif enddo endif wt1 = 1. if (tim(3) .ne. tim(1)) wt1 = (tim(3)-tim(2))/(tim(3)-tim(1)) wt1 = max(0., min(1., wt1)) wt3 = 1. - wt1 dat(locali,localj,2) = dat(locali,localj,1)*wt1 & + dat(locali,localj,3)*wt3 endif write(20,*) "tim(2): ", tim(2) sstreturn = dat(locali,localj,2) write(20,*) "sstreturn; ", sstreturn return end subroutine getvara_verbose (name, ncid, ln, is, ic, dout, s, o) !======================================================================= ! read data ! input: ! name = name of variable to be written ! ncid = iou unit ! ln = length of data ! is = starting point for read in each dimension ! ic = count (or length) for read in each dimension ! s = data scalar ! o = data offset ! output: ! dout = data (default real) !======================================================================= implicit none include 'netcdf.inc' character(*), intent(in) :: name integer, intent(in) :: is(10), ic(10), ln, ncid integer i, iv, nd real, intent(in) :: o, s real dout(ln) real(kind=8) din(ln), offset, scale i = nf_inq_varid (ncid, name, iv) write(20,*) 'ncid: ', ncid write(20,*) 'name: ', name write(20,*) 'var_id: ', iv if (i .ne. nf_noerr) then print*, '==> Warning: netcdf variable ',trim(name), ' not found' return endif scale = 1.0 offset = 0.0 i = nf_inq_varndims(ncid, iv, nd) write(20,*) 'num dims: ', nd call checkerror (i,'getvara nf_inq_varndims '//name) i = nf_get_att_double (ncid, iv, 'add_offset', offset) write(20,*) 'offset: ', offset i = nf_get_att_double (ncid, iv, 'scale_factor', scale) write(20,*) 'scale factor: ',scale i = nf_get_vara_double (ncid, iv, is(1:nd), ic(1:nd), din) call checkerror(i,'getvara '//name) dout(1:ln) = (din(1:ln)*scale + offset)*s + o write(20,*) 'leaving getvara_verbose ..' return end