program downscale_tau
! Purpose
! -------
! Interpolates climatological surface stresses from the Mk3L atmosphere model
! grid on to the grid used by the new, high-resolution ocean model. This
! interpolation is performed in exactly the same fashion as within the coupled
! model.
!
! Usage
! -----
! ./downsale_tau <in_file> <out_file>
!
! where:
!
! in_file file containing climatological surface stresses from the AGCM
! out_file output file
!
! History
! -------
! 2008 Feb 22 Steven Phipps Original version
implicit none
include 'netcdf.inc'
! Define parameters
integer, parameter :: imt = 130, &
jmt = 114, &
nx = 64, &
ny = 56
real, parameter :: pi = 3.14159265, &
radian = 180.0/pi, &
radius = 6.37122e8, &
swldeg = -89.1331111924962
real, dimension(58), parameter :: alatu = (/ -1.57079632342145, &
-1.51519509774365, &
-1.45959388227021, &
-1.40399269429295, &
-1.34839154842973, &
-1.29279046039513, &
-1.23718944745986, &
-1.18158852900217, &
-1.12598772719084, &
-1.07038706785422, &
-1.01478658161207, &
-0.959186305381917, &
-0.903586284425691, &
-0.847986575189126, &
-0.792387249328565, &
-0.736788399559000, &
-0.681190148374214, &
-0.625592661443904, &
-0.569996168916890, &
-0.514401000689571, &
-0.458807647669979, &
-0.403216874610710, &
-0.347629943659605, &
-0.292049100964518, &
-0.236478778882343, &
-0.180929158247486, &
-0.125430272644584, &
-7.013097509103954E-002, &
0.000000000000000E+000, &
7.013097509103954E-002, &
0.125430272644584, &
0.180929158247486, &
0.236478778882343, &
0.292049100964518, &
0.347629943659605, &
0.403216874610710, &
0.458807647669979, &
0.514401000689571, &
0.569996168916890, &
0.625592661443904, &
0.681190148374214, &
0.736788399559000, &
0.792387249328565, &
0.847986575189126, &
0.903586284425691, &
0.959186305381917, &
1.01478658161207, &
1.07038706785422, &
1.12598772719084, &
1.18158852900217, &
1.23718944745986, &
1.29279046039513, &
1.34839154842973, &
1.40399269429295, &
1.45959388227021, &
1.51519509774365, &
1.57079632342145, &
1.62639754909925 /)
! Declare local variables
character(len=80) :: in_file, out_file
integer :: fileid, i, ia, iap1, j, ja, jap1, latdid, latvid, londid, &
lonvid, mondid, month, monvid, status, tauxid, tauyid
real :: sumdy
real(kind=4), dimension(imt-2) :: lonuv
real(kind=4), dimension(jmt-2) :: latuv
real, dimension(jmt) :: dyt, phi
real(kind=4), dimension(nx, ny, 12) :: staux, stauy
real(kind=4), dimension(imt-2, jmt-2, 12) :: otx, oty
! ---------------------
! BEGIN EXECUTABLE CODE
! ---------------------
! ------------------------------------------------------------------------
! The following section calculates the values of PHI, the latitude of each
! velocity gridbox, in exactly the same way as in ocean.f
! ------------------------------------------------------------------------
! Set the meridional dimension of each tracer gridbox in degrees
dyt( 2) = 1.5717758068763
dyt( 3) = 1.5717758068763
dyt( 4) = 1.5873125395736
dyt( 5) = 1.5873125395736
dyt( 6) = 1.5904318409297
dyt( 7) = 1.5904318409297
dyt( 8) = 1.5915276112852
dyt( 9) = 1.5915276112852
dyt( 10) = 1.5920305756336
dyt( 11) = 1.5920305756336
dyt( 12) = 1.5923010877064
dyt( 13) = 1.5923010877064
dyt( 14) = 1.5924626726494
dyt( 15) = 1.5924626726494
dyt( 16) = 1.5925666534402
dyt( 17) = 1.5925666534402
dyt( 18) = 1.5926373800807
dyt( 19) = 1.5926373800807
dyt( 20) = 1.5926875876849
dyt( 21) = 1.5926875876849
dyt( 22) = 1.5927244533195
dyt( 23) = 1.5927244533195
dyt( 24) = 1.5927522680130
dyt( 25) = 1.5927522680130
dyt( 26) = 1.5927737232154
dyt( 27) = 1.5927737232154
dyt( 28) = 1.5927905752292
dyt( 29) = 1.5927905752292
dyt( 30) = 1.5928040087625
dyt( 31) = 1.5928040087625
dyt( 32) = 1.5928148456915
dyt( 33) = 1.5928148456915
dyt( 34) = 1.5928236699541
dyt( 35) = 1.5928236699541
dyt( 36) = 1.5928309049419
dyt( 37) = 1.5928309049419
dyt( 38) = 1.5928368629406
dyt( 39) = 1.5928368629406
dyt( 40) = 1.5928417775583
dyt( 41) = 1.5928417775583
dyt( 42) = 1.5928458254894
dyt( 43) = 1.5928458254894
dyt( 44) = 1.5928491414079
dyt( 45) = 1.5928491414079
dyt( 46) = 1.5928518283288
dyt( 47) = 1.5928518283288
dyt( 48) = 1.5928539649014
dyt( 49) = 1.5928539649014
dyt( 50) = 1.5928556105851
dyt( 51) = 1.5928556105851
dyt( 52) = 1.5928568093168
dyt( 53) = 1.5928568093168
dyt( 54) = 1.5928575920821
dyt( 55) = 1.5928575920821
dyt( 56) = 1.5928579786506
dyt( 57) = 1.5928579786506
dyt( 58) = 1.5928579786506
dyt( 59) = 1.5928579786506
dyt( 60) = 1.5928575920821
dyt( 61) = 1.5928575920821
dyt( 62) = 1.5928568093168
dyt( 63) = 1.5928568093168
dyt( 64) = 1.5928556105851
dyt( 65) = 1.5928556105851
dyt( 66) = 1.5928539649014
dyt( 67) = 1.5928539649014
dyt( 68) = 1.5928518283288
dyt( 69) = 1.5928518283288
dyt( 70) = 1.5928491414079
dyt( 71) = 1.5928491414079
dyt( 72) = 1.5928458254894
dyt( 73) = 1.5928458254894
dyt( 74) = 1.5928417775583
dyt( 75) = 1.5928417775583
dyt( 76) = 1.5928368629406
dyt( 77) = 1.5928368629406
dyt( 78) = 1.5928309049419
dyt( 79) = 1.5928309049419
dyt( 80) = 1.5928236699541
dyt( 81) = 1.5928236699541
dyt( 82) = 1.5928148456915
dyt( 83) = 1.5928148456915
dyt( 84) = 1.5928040087625
dyt( 85) = 1.5928040087625
dyt( 86) = 1.5927905752292
dyt( 87) = 1.5927905752292
dyt( 88) = 1.5927737232154
dyt( 89) = 1.5927737232154
dyt( 90) = 1.5927522680130
dyt( 91) = 1.5927522680130
dyt( 92) = 1.5927244533195
dyt( 93) = 1.5927244533195
dyt( 94) = 1.5926875876849
dyt( 95) = 1.5926875876849
dyt( 96) = 1.5926373800807
dyt( 97) = 1.5926373800807
dyt( 98) = 1.5925666534402
dyt( 99) = 1.5925666534402
dyt(100) = 1.5924626726494
dyt(101) = 1.5924626726494
dyt(102) = 1.5923010877064
dyt(103) = 1.5923010877064
dyt(104) = 1.5920305756336
dyt(105) = 1.5920305756336
dyt(106) = 1.5915276112852
dyt(107) = 1.5915276112852
dyt(108) = 1.5904318409297
dyt(109) = 1.5904318409297
dyt(110) = 1.5873125395736
dyt(111) = 1.5873125395736
dyt(112) = 1.5717758068763
dyt(113) = 1.5717758068763
! Convert DYT to centimetres
dyt = dyt * radius / radian
! Calculate the latitude of each velocity gridbox
phi(1) = swldeg / radian
sumdy = phi(1)
do j = 1, jmt
if (j /= jmt) then
sumdy = sumdy + dyt(j+1) / radius
phi(j+1) = sumdy
end if
end do
! -----------------------------------------------------
! Get the climatological surface stresses from the AGCM
! -----------------------------------------------------
! Get the name of the input file
call getarg(1, in_file)
! Get the climatological surface fluxes from the AGCM
status = nf_open(in_file, nf_nowrite, fileid)
status = nf_inq_varid(fileid, "taux", tauxid)
status = nf_inq_varid(fileid, "tauy", tauyid)
status = nf_get_var_real(fileid, tauxid, staux)
status = nf_get_var_real(fileid, tauyid, stauy)
status = nf_close(fileid)
! ------------------------------------------------------------------------
! The following section of code interpolates the surface stresses from the
! AGCM grid to the new, high-resolution OGCM grid.
!
! These calculations are performed in exactly the same way as in ocntau.f.
! ------------------------------------------------------------------------
! Loop over months
do month = 1, 12
do j = 1, jmt-2
ja = j / 2
jap1 = (j + 1) / 2
do i = 1, imt-2
ia = i / 2
iap1 = (i + 1) / 2
if (i .eq. 1) ia = (imt - 1)/2
if (mod(i, 2) .eq. 0) then
if (mod(j, 2) .eq. 0) then
! i even, j even
otx(i, j, month) = staux(ia, ja, month)
oty(i, j, month) = stauy(ia, ja, month)
else
! i even, j odd
if (j .eq. 1) then
otx(i, j, month) = staux(ia, jap1, month)
oty(i, j, month) = stauy(ia, jap1, month)
else
otx(i, j, month) = &
(cos(alatu(ja+1)) * staux(ia, ja, month) + &
cos(alatu(ja+2)) * staux(ia, jap1, month)) / &
(cos(alatu(ja+1)) + cos(alatu(ja+2)))
oty(i, j, month) = &
(cos(alatu(ja+1)) * stauy(ia, ja, month) + &
cos(alatu(ja+2)) * stauy(ia, jap1, month)) / &
(cos(alatu(ja+1)) + cos(alatu(ja+2)))
end if
endif
else
if (mod(j, 2) .eq. 0) then
! i odd, j even
otx(i, j, month) = 0.5 * (staux(ia, ja, month) + &
staux(iap1, ja, month))
oty(i, j, month) = 0.5 * (stauy(ia, ja, month) + &
stauy(iap1, ja, month))
else
! i odd, j odd
if (j .eq. 1) then
otx(i, j, month) = 0.5 * (staux(ia, jap1, month) + &
staux(iap1, jap1, month))
oty(i, j, month) = 0.5 * (stauy(ia, jap1, month) + &
stauy(iap1, jap1, month))
else
otx(i, j, month) = &
(cos(alatu(ja+1)) * staux(ia, ja, month) + &
cos(alatu(ja+1)) * staux(iap1, ja, month) + &
cos(alatu(ja+2)) * staux(ia, jap1, month) + &
cos(alatu(ja+2)) * staux(iap1, jap1, month)) / &
(2.0 * (cos(alatu(ja+1)) + cos(alatu(ja+2))))
oty(i, j, month) = &
(cos(alatu(ja+1)) * stauy(ia, ja, month) + &
cos(alatu(ja+1)) * stauy(iap1, ja, month) + &
cos(alatu(ja+2)) * stauy(ia, jap1, month) + &
cos(alatu(ja+2)) * stauy(iap1, jap1, month)) / &
(2.0 * (cos(alatu(ja+1)) + cos(alatu(ja+2))))
end if
end if
end if
end do
end do
end do
! --------------------------------------------------------
! Derive the axis data for the interpolated surface fluxes
! --------------------------------------------------------
do i = 1, imt-2
lonuv(i) = (real(i) - 1.0) * 360.0 / float(imt-2)
end do
do j = 1, jmt-2
latuv(j) = phi(j+1) * radian
end do
! --------------------------------
! Save the data to the output file
! --------------------------------
! Get the name of the output file
call getarg(2, out_file)
! (1) Create netCDF file
status = nf_create(out_file, nf_noclobber, fileid)
! (2) Define dimensions
status = nf_def_dim(fileid, "lonuv", imt-2, londid)
status = nf_def_dim(fileid, "latuv", jmt-2, latdid)
status = nf_def_dim(fileid, "month", 12, mondid)
! (3) Define variables
status = nf_def_var(fileid, "lonuv", nf_float, 1, londid, lonvid)
status = nf_def_var(fileid, "latuv", nf_float, 1, latdid, latvid)
status = nf_def_var(fileid, "month", nf_int, 1, mondid, monvid)
status = nf_def_var(fileid, "taux", nf_float, 3, &
(/ londid, latdid, mondid /), tauxid)
status = nf_def_var(fileid, "tauy", nf_float, 3, &
(/ londid, latdid, mondid /), tauyid)
! (4) Put attributes
status = nf_put_att_text(fileid, lonvid, "long_name", 26, &
"Longitude of UV gridpoints")
status = nf_put_att_text(fileid, latvid, "long_name", 25, &
"Latitude of UV gridpoints")
status = nf_put_att_text(fileid, tauxid, "long_name", 17, &
"Zonal wind stress")
status = nf_put_att_text(fileid, tauyid, "long_name", 22, &
"Meridional wind stress")
status = nf_put_att_text(fileid, lonvid, "units", 12, "degrees_east")
status = nf_put_att_text(fileid, latvid, "units", 13, "degrees_north")
status = nf_put_att_text(fileid, tauxid, "units", 10, "dynes/cm^2")
status = nf_put_att_text(fileid, tauyid, "units", 10, "dynes/cm^2")
status = nf_put_att_text(fileid, monvid, "units", 6, "months")
status = nf_put_att_text(fileid, lonvid, "point_spacing", 4, "even")
status = nf_put_att_text(fileid, latvid, "point_spacing", 6, "uneven")
status = nf_put_att_text(fileid, monvid, "point_spacing", 4, "even")
status = nf_put_att_text(fileid, lonvid, "modulo", 1, " ")
status = nf_put_att_text(fileid, monvid, "modulo", 1, " ")
! (5) Exit define mode
status = nf_enddef(fileid)
! (6) Write data
status = nf_put_var_real(fileid, lonvid, lonuv)
status = nf_put_var_real(fileid, latvid, latuv)
status = nf_put_var_int(fileid, monvid, (/ 1, 2, 3, 4, 5, 6, &
7, 8, 9, 10, 11, 12 /))
status = nf_put_var_real(fileid, tauxid, otx)
status = nf_put_var_real(fileid, tauyid, oty)
! (7) Close output file
status = nf_close(fileid)
end program downscale_tau