! Generated by TAPENADE (INRIA, Tropics team)
! Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!
! $Id: ice_flux_in.F,v 1.14 2004/03/01 16:56:23 eclare Exp $
!=======================================================================
!
!BOP
!
! !MODULE: ice_flux_in - reads and interpolates input forcing data
!
! !DESCRIPTION:
!
! Reads and interpolates forcing data for atmospheric and oceanic quantities.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !INTERFACE:
!
MODULE
ICE_FLUX_IN
USE
ice_calendarUSE
ice_constantsUSE
ice_diagnosticsUSE
ice_domainUSE
ice_fileunitsUSE
ice_fluxUSE
ice_gridUSE
ice_historyUSE
ice_kinds_modUSE
ice_model_sizeUSE
ice_oceanUSE
ice_read_writeUSE
ice_workCHARACTER
(len
=char_len_long) ::atm_data_dir
,ocn_data_dir
REAL(KIND=DBL_KIND)
:: cldf(ilo:ihi, jlo:jhi)REAL(KIND=DBL_KIND)
::c1intp
,c2intp
,ftime
INTEGER(KIND=INT_KIND)
::fyear
INTEGER(KIND=INT_KIND)
::oldrecnum
=0
CHARACTER
(len
=char_len_long) ::flw_file
,fsw_file
,height_file
,&
&
humid_file
,pott_file
,rain_file
,rhoa_file
,sss_file
,sst_file
,&
&
tair_file
,uwind_file
,vwind_file
REAL(KIND=DBL_KIND)
:: cldf_data(ilo:ihi, jlo:jhi,2
), flw_data(ilo:&
&
ihi, jlo:jhi,2
), fsnow_data(ilo:ihi, jlo:jhi,2
), fsw_data(ilo:ihi,&
&
jlo:jhi,2
), pott_data(ilo:ihi, jlo:jhi,2
), qa_data(ilo:ihi, jlo:jhi&
&
,2
), rhoa_data(ilo:ihi, jlo:jhi,2
), sss_data(ilo:ihi, jlo:jhi,2
),&
&
sst_data(ilo:ihi, jlo:jhi,2
), tair_data(ilo:ihi, jlo:jhi,2
),&
&
uatm_data(ilo:ihi, jlo:jhi,2
), vatm_data(ilo:ihi, jlo:jhi,2
),&
&
zlvl_data(ilo:ihi, jlo:jhi,2
)INTEGER(KIND=INT_KIND)
::fyear_final
,fyear_init
,ycycle
CONTAINS
SUBROUTINE
INIT_GETFLUX
()IMPLICIT NONE
INTEGER
::arg1
INTRINSIC
MOD
!
! !DESCRIPTION:
!
! Determines the current and final years of the forcing cycle based on
! namelist input, initializes the forcing data filenames, and initializes
! surface temperature and salinity from data. \\
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! last year in forcing cycle
fyear_final = fyear_init + ycycle -1
! current year
arg1 = nyr - year_init fyear = fyear_init +MOD
(arg1, ycycle)!
IF
(my_task .EQ. master_task)THEN
WRITE
(nu_diag,*
)' Initial forcing data year = '
, fyear_initWRITE
(nu_diag,*
)' Current forcing data year = '
, fyearWRITE
(nu_diag,*
)' Final forcing data year = '
, fyear_finalEND IF
!
!ech-later call Tair_limit
!
CALL
NCAR_FILES
(fyear)END
SUBROUTINE
INIT_GETFLUX
SUBROUTINE
GETFLUX
()IMPLICIT NONE
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
INTEGER
::arg1
INTRINSIC
MOD
!
! !DESCRIPTION:
!
! Get forcing data and interpolate as necessary!
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! current year
arg1 = nyr - year_init fyear = fyear_init +MOD
(arg1, ycycle)!
! forcing time
ftime = time! for restarting
!-------------------------------------------------------------------
! Read and interpolate annual climatologies of SSS and SST.
! Restore model SST to data.
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Read and interpolate atmospheric data
!-------------------------------------------------------------------
time_forc = ftime!
!
! call sss_sst_restore
!
CALL
NCAR_BULK_DAT
()CALL
PREPARE_FORCING
()END
SUBROUTINE
GETFLUX
SUBROUTINE
NCAR_FILES
(yr
)IMPLICIT NONE
INTEGER(KIND=INT_KIND),INTENT(IN)
::yr
!
! !DESCRIPTION:
!
! This subroutine is based on the LANL file naming conventions.
! Edit for other directory structures or filenames.
! Note: The year number in these filenames does not matter, because
! subroutine file\_year will insert the correct year.
!
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
! current forcing year
!
!EOP
!
fsw_file ='/home/jkim/atm/swdn.1997.dat'
CALL
FILE_YEAR
(fsw_file, yr)!
flw_file ='/home/jkim/atm/cldf.1997.dat'
CALL
FILE_YEAR
(flw_file, yr)!
rain_file ='/home/jkim/atm/prec.1997.dat'
CALL
FILE_YEAR
(rain_file, yr)!
uwind_file ='/home/jkim/atm/u_10.1997.dat'
CALL
FILE_YEAR
(uwind_file, yr)!
vwind_file ='/home/jkim/atm/v_10.1997.dat'
CALL
FILE_YEAR
(vwind_file, yr)!
tair_file ='/home/jkim/atm/t_10.1997.dat'
CALL
FILE_YEAR
(tair_file, yr)!
humid_file ='/home/jkim/atm/q_10.1997.dat'
CALL
FILE_YEAR
(humid_file, yr)!
rhoa_file ='/home/jkim/atm/dn10.1997.dat'
CALL
FILE_YEAR
(rhoa_file, yr)! master_task
!
IF
(my_task .EQ. master_task)THEN
WRITE
(nu_diag,*
)''
WRITE
(nu_diag,*
)'Forcing data year = '
, fyearWRITE
(nu_diag,*
)'Atmospheric data files:'
WRITE
(nu_diag,*
) fsw_fileWRITE
(nu_diag,*
) flw_fileWRITE
(nu_diag,*
) rain_fileWRITE
(nu_diag,*
) uwind_fileWRITE
(nu_diag,*
) vwind_fileWRITE
(nu_diag,*
) tair_fileWRITE
(nu_diag,*
) humid_fileWRITE
(nu_diag,*
) rhoa_fileEND IF
END
SUBROUTINE
NCAR_FILES
SUBROUTINE
NCAR_BULK_DAT
()IMPLICIT NONE
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
INTEGER
::arg1
DOUBLE PRECISION
::arg10
INTEGER(KIND=INT_KIND)
::i
,imx
,ipx
,ixx
,j
,maxrec
,midmonth
,&
&
recnum
,recslot
LOGICAL(KIND=LOG_KIND)
::read6
,readm
REAL(KIND=DBL_KIND)
::sec6hr
INTRINSIC
MOD
,INT
,DBLE
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! record numbers for neighboring months
! record number
! maximum record number
! spline slot for current record
! middle day of month
!
! number of seconds in 6 hours
!
!-------------------------------------------------------------------
! monthly data
!
! Assume that monthly data values are located in the middle of the
! month.
!-------------------------------------------------------------------
!
! data is given on 15th of every month
! Compute record numbers for surrounding months
midmonth =15
! midmonth = fix(p5 * dble(daymo(month))) ! exact middle
!
maxrec =12
arg1 = month + maxrec -2
imx =MOD
(arg1, maxrec) +1
ipx =MOD
(month, maxrec) +1
! other two points will be used
IF
(mday .GE. midmonth) imx =99
! Determine whether interpolation will use values 1:2 or 2:3
! recslot = 2 means we use values 1:2, with the current value (2)
! in the second slot
! recslot = 1 means we use values 2:3, with the current value (2)
! in the first slot
IF
(mday .LT. midmonth) ipx =99
!
! latter half of month
recslot =1
! first half of month
! Find interpolation coefficients
IF
(mday .LT. midmonth) recslot =2
!
CALL
INTERP_COEFF_MONTHLY
(recslot)! Read 2 monthly values
!
readm =.false.
IF
(istep .EQ.1
.OR. mday .EQ. midmonth .AND. sec .EQ.0
) readm =&
&
.true.
!
CALL
READ_DATA
(readm,0
, fyear, imx, month, ipx, maxrec, fsw_file,&
&
fsw_data)CALL
READ_DATA
(readm,0
, fyear, imx, month, ipx, maxrec, flw_file,&
&
cldf_data)CALL
READ_DATA
(readm,0
, fyear, imx, month, ipx, maxrec, rain_file,&
&
fsnow_data)!
CALL
INTERPOLATE_DATA
(fsw_data, fsw)CALL
INTERPOLATE_DATA
(cldf_data, cldf)CALL
INTERPOLATE_DATA
(fsnow_data, fsnow)!-------------------------------------------------------------------
! 6-hourly data
!
! Assume that the 6-hourly value is located at the end of the
! 6-hour period. This is the convention for NCEP reanalysis data.
! E.g. record 1 gives conditions at 6 am GMT on 1 January.
!-------------------------------------------------------------------
!
!
! seconds in 6 hours
sec6hr = secday/c4! 365*4
! current record number
maxrec =1460
!
! Compute record numbers for surrounding data (2 on each side)
arg10 =DBLE
(sec)/sec6hr recnum =4
*INT
(yday) -3
+INT
(arg10)!
!
arg1 = recnum + maxrec -2
imx =MOD
(arg1, maxrec) +1
! Compute interpolation coefficients
! If data is located at the end of the time interval, then the
! data value for the current record goes in slot 2
arg1 = recnum -1
ixx =MOD
(arg1, maxrec) +1
! ipx = mod(recnum, maxrec) + 1
!
!
recslot =2
ipx =99
CALL
INTERP_COEFF
(recnum, recslot, sec6hr)! Read
!
read6 =.false.
IF
(istep .EQ.1
.OR. oldrecnum .NE. recnum) read6 =.true.
!
CALL
READ_DATA
(read6,0
, fyear, imx, ixx, ipx, maxrec, tair_file,&
&
tair_data)CALL
READ_DATA
(read6,0
, fyear, imx, ixx, ipx, maxrec, uwind_file,&
&
uatm_data)CALL
READ_DATA
(read6,0
, fyear, imx, ixx, ipx, maxrec, vwind_file,&
&
vatm_data)CALL
READ_DATA
(read6,0
, fyear, imx, ixx, ipx, maxrec, rhoa_file,&
&
rhoa_data)CALL
READ_DATA
(read6,0
, fyear, imx, ixx, ipx, maxrec, humid_file,&
&
qa_data)! Interpolate
!
CALL
INTERPOLATE_DATA
(tair_data, tair)CALL
INTERPOLATE_DATA
(uatm_data, uatm)CALL
INTERPOLATE_DATA
(vatm_data, vatm)CALL
INTERPOLATE_DATA
(rhoa_data, rhoa)CALL
INTERPOLATE_DATA
(qa_data, qa)! Save record number
!
oldrecnum = recnumEND
SUBROUTINE
NCAR_BULK_DAT
SUBROUTINE
READ_DATA
(flag
,recd
,yr
,imx
,ixx
,ipx
,maxrec
,data_file
&
&
,field_data
)IMPLICIT NONE
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
CHARACTER
(len
=char_len_long) ::data_file
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi,2
),INTENT(OUT)
::&
&
field_data
LOGICAL(KIND=LOG_KIND),INTENT(IN)
::flag
INTEGER(KIND=INT_KIND),INTENT(IN)
::imx
INTEGER(KIND=INT_KIND),INTENT(IN)
::ipx
INTEGER(KIND=INT_KIND),INTENT(IN)
::ixx
INTEGER(KIND=INT_KIND),INTENT(IN)
::maxrec
INTEGER(KIND=INT_KIND),INTENT(IN)
::recd
INTEGER(KIND=INT_KIND),INTENT(IN)
::yr
INTEGER(KIND=INT_KIND)
::arg1
REAL(KIND=DBL_KIND)
:: field_dummy(ilo:ihi, jlo:jhi)INTEGER(KIND=INT_KIND)
::arg
,n2
,n4
,nbits
,nrec
LOGICAL(KIND=LOG_KIND)
::diag
,scatter
!
! !DESCRIPTION:
!
! If data is at the beginning of a one-year record, get data from
! the previous year.
! If data is at the end of a one-year record, get data from the
! following year.
! If no earlier data exists (beginning of fyear\_init), then \\
! (1) For monthly data, get data from the end of fyear\_final. \\
! (2) For more frequent data, let the imx value equal the
! first value of the year. \\
! If no later data exists (end of fyear\_final), then \\
! (1) For monthly data, get data from the beginning of fyear\_init. \\
! (2) For more frequent data, let the ipp and, if necessary, ipx
! values equal the last value of the year. \\
! In other words, we assume persistence when daily or 6-hourly
! data is missing, and we assume periodicity when monthly data
! is missing. \\
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !
!
! !INPUT/OUTPUT PARAMETERS:
!
!
! baseline record number
! year of forcing data
! record numbers of 3 data values
! relative to recd
! maximum record value
!
! 2 values needed for interpolation
!
!EOP
!
! data file to be read
!
! = 32 for single precision, 64 for double
! record number to read
! like imx and ipx, but
! adjusted at beginning and end of data
! value of time argument in field_data
!
!
! double precision data
nbits =64
!
! scatter data to processors
scatter =.true.
! write diagnostic info
diag =.false.
!
!! debugging
IF
(istep1 .GT. check_step) diag =.true.
! flag
!
! if (my_task==master_task .and. (diag)) then
! write(nu_diag,*) ' ',data_file
! endif
!
IF
(flag)THEN
!-----------------------------------------------------------------
! Initialize record counters
! (n2, n4 will change only at the very beginning or end of
! a forcing cycle.)
!-----------------------------------------------------------------
n2 = imx n4 = ipx!-----------------------------------------------------------------
! read data
!-----------------------------------------------------------------
arg =0
! imx ne 99
! always read ixx data from data file for current year
!
IF
(imx .NE.99
)THEN
! ixx <= 1
! currently in first half of data interval
IF
(ixx .LE.1
)THEN
! yr > fyear_init
IF
(yr .GT. fyear_init)THEN
! get data from previous year
arg1 = yr -1
CALL
FILE_YEAR
(data_file, arg1)ELSE
IF
(maxrec .GT.12
)THEN
! yr = fyear_init, no prior data exists
! extrapolate from first record
IF
(ixx .EQ.1
) n2 = ixxELSE
! go to end of fyear_final
CALL
FILE_YEAR
(data_file, fyear_final)END IF
END IF
CALL
ICE_OPEN
(nu_forcing, data_file, nbits)!
arg =1
nrec = recd + n2!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter,&
&
diag) field_data(:, :, arg) = field_dummy(:, :)!
IF
(ixx .EQ.1
.AND. my_task .EQ. master_task)CLOSE
(nu_forcing)&
&
END IF
!
CALL
FILE_YEAR
(data_file, yr)CALL
ICE_OPEN
(nu_forcing, data_file, nbits)!
arg = arg +1
nrec = recd + ixx!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter, diag&
&
) field_data(:, :, arg) = field_dummy(:, :)! ipx ne 99
!
IF
(ipx .NE.99
)THEN
! ixx = maxrec
! currently in latter half of data interval
IF
(ixx .EQ. maxrec)THEN
! yr < fyear_final
IF
(yr .LT. fyear_final)THEN
! get data from following year
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing) arg1 = yr +1
CALL
FILE_YEAR
(data_file, arg1)CALL
ICE_OPEN
(nu_forcing, data_file, nbits)ELSE
IF
(maxrec .GT.12
)THEN
! yr = fyear_final, no more data exists
! extrapolate from ixx
n4 = ixxELSE
! go to beginning of fyear_init
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing)CALL
FILE_YEAR
(data_file, fyear_init)CALL
ICE_OPEN
(nu_forcing, data_file, nbits)END IF
END IF
!
arg = arg +1
nrec = recd + n4!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter,&
&
diag) field_data(:, :, arg) = field_dummy(:, :)END IF
!
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing)END IF
END
SUBROUTINE
READ_DATA
SUBROUTINE
FILE_YEAR
(data_file
,yr
)IMPLICIT NONE
CHARACTER
(len
=char_len_long),INTENT(INOUT)
::data_file
INTEGER(KIND=INT_KIND),INTENT(IN)
::yr
INTEGER(KIND=INT_KIND)
::i
CHARACTER
(len
=char_len_long) ::tmpname
INTRINSIC
INDEX
!
! !DESCRIPTION:
!
! Construct the correct name of the atmospheric data file
! to be read, given the year and assuming the naming convention
! that filenames end with 'yyyy.dat'.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!
!EOP
!
!
!
!
i =INDEX
(data_file,'.dat'
) -5
tmpname = data_fileWRITE
(data_file,'(a,i4.4,a)'
) tmpname(1
:i), yr,'.dat'
END
SUBROUTINE
FILE_YEAR
SUBROUTINE
INTERP_COEFF
(recnum
,recslot
,secint
)IMPLICIT NONE
INTEGER(KIND=INT_KIND),INTENT(IN)
::recnum
INTEGER(KIND=INT_KIND),INTENT(IN)
::recslot
REAL(KIND=DBL_KIND),INTENT(IN)
::secint
REAL(KIND=DBL_KIND)
,PARAMETER
::secyr
=c365*secdayREAL(KIND=DBL_KIND)
::arg1
REAL(KIND=DBL_KIND)
::t1
,t2
,tt
INTRINSIC
ABS
,INT
,DBLE
!
! !DESCRIPTION:
!
! Compute coefficients for interpolating data to current time step.
! Works for any data interval that divides evenly into a 365-day
! year (daily, 6-hourly, etc.)
! Use interp\_coef\_monthly for monthly data.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !USES:
!
! record number for current data value
! spline slot for current record
!
! seconds in data interval
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! seconds in a 365-day year
!
! seconds elapsed in current year
! seconds elapsed at data points
!
!
!jkim tt = mod(ftime,secyr)
! Find neighboring times
arg1 = ftime/secyr tt = ftime -INT
(arg1)*secyr! Compute coefficients
IF
(recslot .EQ.2
)THEN
! current record goes in slot 2 (NCEP)
t2 =DBLE
(recnum)*secint! - 1 interval
t1 = t2 - secintELSE
! recslot = 1
t1 =DBLE
(recnum)*secint! + 1 interval
t2 = t1 + secintEND IF
IF
((t2-tt)/(t2-t1) .GE.0.
)THEN
c1intp = (t2-tt)/(t2-t1)ELSE
c1intp = -((t2-tt)/(t2-t1))END IF
c2intp = c1 - c1intpEND
SUBROUTINE
INTERP_COEFF
SUBROUTINE
PREPARE_FORCING
()IMPLICIT NONE
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
REAL(KIND=DBL_KIND)
::arg1
INTEGER(KIND=INT_KIND)
::i
,j
REAL(KIND=DBL_KIND)
::workx
,worky
REAL
::x1
INTRINSIC
COS
,EXP
,SIN
,MAX
,MIN
,SQRT
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
DO
j=jlo,jhi!-----------------------------------------------------------------
! Fix interpolation errors
!-----------------------------------------------------------------
DO
i=ilo,ihiIF
(fsw(i, j) .LT. c0)THEN
fsw(i, j) = c0ELSE
fsw(i, j) = fsw(i, j)END IF
IF
(cldf(i, j) .GT. c1)THEN
x1 = c1ELSE
x1 = cldf(i, j)END IF
IF
(x1 .LT. c0)THEN
cldf(i, j) = c0ELSE
cldf(i, j) = x1END IF
IF
(fsnow(i, j) .LT. c0)THEN
fsnow(i, j) = c0ELSE
fsnow(i, j) = fsnow(i, j)END IF
IF
(rhoa(i, j) .LT. c0)THEN
rhoa(i, j) = c0ELSE
rhoa(i, j) = rhoa(i, j)END IF
IF
(qa(i, j) .LT. c0)THEN
qa(i, j) = c0ELSE
qa(i, j) = qa(i, j)END IF
!
! as in the dummy atm (latm)
qa(i, j) = qa(i, j)*0.94_dbl_kind
! as in the dummy atm (latm)
!-----------------------------------------------------------------
! Compute other fields needed by model
!-----------------------------------------------------------------
fsw(i, j) = fsw(i, j)*0.92_dbl_kind
!
!
zlvl(i, j) = c10! wind speed, m/s
arg1 = uatm(i, j)**2
+ vatm(i, j)**2
wind(i, j) =SQRT
(arg1)! divide shortwave into spectral bands
pott(i, j) = tair(i, j)!
! visible direct
swvdr(i, j) = fsw(i, j)*.28
! visible diffuse
swvdf(i, j) = fsw(i, j)*.24
! near IR direct
swidr(i, j) = fsw(i, j)*.31
! near IR diffuse
! as in the dummy atm (latm)
! longwave as in Parkinson and Washington (1979)
swidf(i, j) = fsw(i, j)*.17
!
!!! downward longwave !!!
! determine whether precip is rain or snow
arg1 = -(7.77e-4_dbl_kind
*(tffresh-tair(i, j))**2
) flw(i, j) = stefan_boltzmann*tair(i, j)**4
*(c1-0.261_dbl_kind
*&
&
EXP
(arg1))*(c1+0.275_dbl_kind
*cldf(i, j))!
! mm/month -> kg/m^2 s
fsnow(i, j) = fsnow(i, j)/2.592e+06_dbl_kind
frain(i, j) = c0!-----------------------------------------------------------------
! rotate zonal/meridional vectors to local coordinates
! Vector fields come in on T grid, but are oriented geographically
! need to rotate to pop-grid FIRST using ANGLET
! then interpolate to the U-cell centers (otherwise we
! interpolate across the pole)
! use ANGLET which is on the T grid !
! atmo variables are needed in T cell centers in subroutine stability,
! and are interpolated to the U grid later as necessary
!-----------------------------------------------------------------
! if (Tair(i,j) >= Tffresh) then
IF
(tair(i, j) .GE. tffresh)THEN
frain(i, j) = fsnow(i, j) fsnow(i, j) = c0END IF
!
! wind velocity, m/s
workx = uatm(i, j) worky = vatm(i, j)! convert to POP grid
! note uatm, vatm, wind
uatm(i, j) = workx*COS
(anglet(i, j)) + worky*SIN
(anglet(i, j))! are on the T-grid here
vatm(i, j) = worky*COS
(anglet(i, j)) - workx*SIN
(anglet(i, j))END DO
END DO
END
SUBROUTINE
PREPARE_FORCING
SUBROUTINE
SSS_SST_RESTORE
()IMPLICIT NONE
INTEGER
::arg1
INTEGER(KIND=INT_KIND)
::i
,imx
,ipx
,j
,maxrec
,midmonth
,recslot
LOGICAL(KIND=LOG_KIND)
::readm
REAL(KIND=DBL_KIND)
:: sstdat(ilo:ihi, jlo:jhi)REAL(KIND=DBL_KIND)
::trestore
INTRINSIC
MOD
!
! !DESCRIPTION:
!
! Interpolate monthly sss, sst data to timestep.
! Restore sst computed by ice model to data.
!
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! record numbers for neighboring months
! maximum record number
! spline slot for current record
! middle day of month
!
! data value toward which SST is restored
!
!
!
!jkim sss_file = trim(ocn_data_dir)//'sss_Lev.mm'
!jkim sst_file = trim(ocn_data_dir)//'sst_Lev.mm'
!
! Trestore = dt ! STRONG restoring
! Trestore = 1.296e6 ! 15 days
! Trestore = 2.592e6 ! 30 days
! 90 days
trestore =7.776e6
!-------------------------------------------------------------------
! monthly data
!
! Assume that monthly data values are located in the middle of the
! month.
!-------------------------------------------------------------------
!
IF
(my_task .EQ. master_task .AND. istep .EQ.1
)THEN
WRITE
(nu_diag,*
)''
WRITE
(nu_diag,*
)'SSS and SST data interpolated to timestep:'
WRITE
(nu_diag,*
) sss_fileWRITE
(nu_diag,*
) sst_fileWRITE
(nu_diag,*
)'SST restoring timescale = '
, trestore/secday,&
&
' days'
END IF
!
! data is given on 15th of every month
! Compute record numbers for surrounding months
midmonth =15
! midmonth = fix(p5 * dble(daymo(month))) ! exact middle
!
maxrec =12
arg1 = month + maxrec -2
imx =MOD
(arg1, maxrec) +1
ipx =MOD
(month, maxrec) +1
! other two points will be used
IF
(mday .GE. midmonth) imx =99
! Determine whether interpolation will use values 1:2 or 2:3
! recslot = 2 means we use values 1:2, with the current value (2)
! in the second slot
! recslot = 1 means we use values 2:3, with the current value (2)
! in the first slot
IF
(mday .LT. midmonth) ipx =99
!
! latter half of month
recslot =1
! first half of month
! Find interpolation coefficients
IF
(mday .LT. midmonth) recslot =2
!
CALL
INTERP_COEFF_MONTHLY
(recslot)! Read 2 monthly values
!
readm =.false.
IF
(istep .EQ.1
.OR. mday .EQ. midmonth .AND. sec .EQ.0
) readm =&
&
.true.
!
CALL
READ_CLIM_DATA
(readm,0
, imx, month, ipx, sss_file, sss_data)CALL
READ_CLIM_DATA
(readm,0
, imx, month, ipx, sst_file, sst_data)! Interpolate
!
CALL
INTERPOLATE_DATA
(sss_data, sss)CALL
INTERPOLATE_DATA
(sst_data, sstdat)!-------------------------------------------------------------------
! Restore sst to data
!-------------------------------------------------------------------
!
DO
j=jlo,jhiDO
i=ilo,ihiIF
(sss(i, j) .LT. c0) sss(i, j) = c0 sst(i, j) = sst(i, j) + (sstdat(i, j)-sst(i, j))*dt/trestoreEND DO
END DO
!
CALL
COMPLETE_GETFLUX_OCN
()END
SUBROUTINE
SSS_SST_RESTORE
SUBROUTINE
INTERP_COEFF_MONTHLY
(recslot
)IMPLICIT NONE
INTEGER(KIND=INT_KIND),INTENT(IN)
::recslot
REAL(KIND=DBL_KIND)
::arg1
REAL(KIND=INT_KIND)
:: daymid(0
:13
)REAL(KIND=DBL_KIND)
::t1
,t2
,tt
INTRINSIC
INT
!
! !DESCRIPTION:
!
! Compute coefficients for interpolating monthly data to current time step.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
! slot (1 or 2) for current record
!
!EOP
!
! seconds elapsed in current year
! seconds elapsed at month midpoint
!
! month mid-points
!
! time frame ends 0 sec into day 15
daymid(1
:13
) =14.
! Dec 15, 0 sec
! make time cyclic
daymid(0
) = -17.
!
!jkim tt = mod(ftime/secday,c365)
! Find neighboring times
arg1 = ftime/secday/c365 tt = ftime/secday -INT
(arg1)*c365! Compute coefficients
IF
(recslot .EQ.2
)THEN
! first half of month
! midpoint, current month
t2 = daycal(month) + daymid(month)IF
(month .EQ.1
)THEN
! Dec 15 (0 sec)
t1 = daymid(0
)ELSE
! midpoint, previous month
t1 = daycal(month-1
) + daymid(month-1
)END IF
ELSE
! second half of month
! midpoint, current month
t1 = daycal(month) + daymid(month)! day 15 of next month (0 sec)
t2 = daycal(month+1
) + daymid(month+1
)END IF
c1intp = (t2-tt)/(t2-t1) c2intp = c1 - c1intpEND
SUBROUTINE
INTERP_COEFF_MONTHLY
SUBROUTINE
READ_CLIM_DATA
(readflag
,recd
,imx
,ixx
,ipx
,data_file
,&
&
field_data
)IMPLICIT NONE
CHARACTER
(len
=char_len_long),INTENT(IN)
::data_file
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi,2
),INTENT(OUT)
::&
&
field_data
INTEGER(KIND=INT_KIND),INTENT(IN)
::imx
INTEGER(KIND=INT_KIND),INTENT(IN)
::ipx
INTEGER(KIND=INT_KIND),INTENT(IN)
::ixx
LOGICAL(KIND=LOG_KIND),INTENT(IN)
::readflag
INTEGER(KIND=INT_KIND),INTENT(IN)
::recd
REAL(KIND=DBL_KIND)
:: field_dummy(ilo:ihi, jlo:jhi)INTEGER(KIND=INT_KIND)
::arg
,nbits
,nrec
LOGICAL(KIND=LOG_KIND)
::diag
,scatter
!
! !DESCRIPTION:
!
! Read data needed for interpolation, as in read\_data.
! Assume a one-year cycle of climatological data, so that there is
! no need to get data from other years or to extrapolate data beyond
! the forcing time period.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
! William H. Lipscomb, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!
! baseline record number
! record numbers of 3 data values
! relative to recd
!
!
! 2 values needed for interpolation
!
!EOP
!
! = 32 for single precision, 64 for double
! record number to read
! value of time argument in field_data
!
!
! double precision data
nbits =64
!
! scatter data to processors
scatter =.true.
! write diagnostic info
diag =.false.
!
!! debugging
IF
(istep1 .GT. check_step) diag =.true.
!
IF
(my_task .EQ. master_task .AND. diag)WRITE
(nu_diag,*
)' '
,&
&
data_file! readflag
!
IF
(readflag)THEN
!-----------------------------------------------------------------
! read data
!-----------------------------------------------------------------
CALL
ICE_OPEN
(nu_forcing, data_file, nbits)!
arg =0
IF
(imx .NE.99
)THEN
arg =1
nrec = recd + imx!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter,&
&
diag) field_data(:, :, arg) = field_dummy(:, :)END IF
!
arg = arg +1
nrec = recd + ixx!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter, diag&
&
) field_data(:, :, arg) = field_dummy(:, :)!
IF
(ipx .NE.99
)THEN
arg = arg +1
nrec = recd + ipx!jkim call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim & 'rda8', scatter, diag)
CALL
ICE_READ
(nu_forcing, nrec, field_dummy,'rda8'
, scatter,&
&
diag) field_data(:, :, arg) = field_dummy(:, :)END IF
!
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing)END IF
END
SUBROUTINE
READ_CLIM_DATA
SUBROUTINE
INTERPOLATE_DATA
(field_data
,field
)IMPLICIT NONE
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::&
&
field
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi,2
),INTENT(IN)
::&
&
field_data
INTEGER(KIND=INT_KIND)
::i
,j
!
! !DESCRIPTION:
!
! Linear interpolation
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
! 2 values used for interpolation
!
! interpolated field
!
!EOP
!
!
DO
j=jlo,jhiDO
i=ilo,ihi field(i, j) = c1intp*field_data(i, j,1
) + c2intp*field_data(i,&
&
j,2
)END DO
END DO
END
SUBROUTINE
INTERPOLATE_DATA
SUBROUTINE
SSS_CLIM
()IMPLICIT NONE
INTEGER(KIND=INT_KIND)
::i
,j
,k
,nbits
LOGICAL(KIND=LOG_KIND)
::diag
,scatter
!
! !DESCRIPTION:
!
! Creates annual mean climatology for Levitus sss from 12-month climatology.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
!
scatter =.true.
diag =.false.
!
! double precision data
nbits =64
!
!jkim sss_file = trim(ocn_data_dir)//'sss_Lev.mm'
!
IF
(my_task .EQ. master_task)THEN
WRITE
(nu_diag,*
)''
WRITE
(nu_diag,*
)'SSS climatology computed from:'
WRITE
(nu_diag,*
) sss_fileEND IF
!
IF
(my_task .EQ. master_task)CALL
ICE_OPEN
(nu_forcing, sss_file,&
&
nbits)!-------------------------------------------------------------------
! create surface salinity climatology from monthly data
!-------------------------------------------------------------------
!
!
DO
j=jlo,jhiDO
i=ilo,ihi worka(i, j) = c0END DO
END DO
!
! loop over 12 months
DO
k=1
,12
!
CALL
ICE_READ
(nu_forcing, k, worka,'rda8'
, scatter, diag)DO
j=jlo,jhiDO
i=ilo,ihi sss(i, j) = worka(i, j) + sss(i, j)END DO
END DO
END DO
!
! k
!
DO
j=jlo,jhiDO
i=ilo,ihi! annual average salinity
sss(i, j) = sss(i, j)/c12IF
(sss(i, j) .LT. c0) sss(i, j) = c0END DO
END DO
! close file
!
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing)!
CALL
COMPLETE_GETFLUX_OCN
()END
SUBROUTINE
SSS_CLIM
SUBROUTINE
COMPLETE_GETFLUX_OCN
()IMPLICIT NONE
INTEGER(KIND=INT_KIND)
::i
,j
INTRINSIC
MAX
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!-------------------------------------------------------------------
! Compute freezing temperature based on salinity
! Make sure SST is not below Tf.
!-------------------------------------------------------------------
!
DO
j=jlo,jhiDO
i=ilo,ihi! freezing temp, deg C
tf(i, j) = -(depresst*sss(i, j))IF
(sst(i, j) .LT. tf(i, j))THEN
sst(i, j) = tf(i, j)ELSE
sst(i, j) = sst(i, j)END IF
END DO
END DO
END
SUBROUTINE
COMPLETE_GETFLUX_OCN
SUBROUTINE
SST_IC
()IMPLICIT NONE
INTEGER(KIND=INT_KIND)
::i
,j
,k
,nbits
LOGICAL(KIND=LOG_KIND)
::diag
,scatter
INTRINSIC
MAX
!
! !DESCRIPTION:
!
! Reads sst data for current month, and adjusts sst based on freezing
! temperature. Does not interpolate.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
!
scatter =.true.
diag =.false.
!
! double precision data
!-----------------------------------------------------------------
! Shea, Trenberth and Reynolds SST
!-----------------------------------------------------------------
nbits =64
!
! sst_file = '/n/filer3/climate/eclare/sst.sss/sst.mm.100x116.da'
!jkim sst_file = trim(ocn_data_dir)//'sst_Lev.mm'
!
IF
(my_task .EQ. master_task)THEN
WRITE
(nu_diag,*
)''
WRITE
(nu_diag,*
)'SST initial condition:'
WRITE
(nu_diag,*
) sst_fileEND IF
!
IF
(my_task .EQ. master_task)CALL
ICE_OPEN
(nu_forcing, sst_file,&
&
nbits)!
CALL
ICE_READ
(nu_forcing, month, sst,'rda8'
, scatter, diag)!
IF
(my_task .EQ. master_task)CLOSE
(nu_forcing)!
DO
j=jlo,jhi! Make sure sst is not less than Tf
DO
i=ilo,ihiIF
(sst(i, j) .LT. tf(i, j))THEN
sst(i, j) = tf(i, j)ELSE
sst(i, j) = sst(i, j)END IF
END DO
END DO
END
SUBROUTINE
SST_IC
END MODULE
ICE_FLUX_IN