program downscale_athf ! Purpose ! ------- ! Interpolates climatological surface fluxes 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. ! ! The fluxes that can be interpolated are the surface heat flux [ATHF] and the ! surface salinity tendency [ATSF]. ! ! Usage ! ----- ! ./downsale_athf ! ! where: ! ! imsl_file file containing the array IMSL ! in_file file containing climatological surface fluxes from the AGCM ! var_name variable containing climatological surface fluxes ! out_file output file ! ! History ! ------- ! 2008 Feb 21 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, & eradsq = 6.37122e6*6.37122e6 real, dimension(28), parameter :: w = (/ 2.323855375777293E-003, & 5.402522246016385E-003, & 8.469063163306595E-003, & 1.150982434038331E-002, & 1.451508927802266E-002, & 1.747551291140060E-002, & 2.038192988240281E-002, & 2.322535156256496E-002, & 2.599698705839140E-002, & 2.868826847382287E-002, & 3.129087674730992E-002, & 3.379676711561259E-002, & 3.619819387231452E-002, & 3.848773425924692E-002, & 4.065831138474480E-002, & 4.270321608466642E-002, & 4.461612765269255E-002, & 4.639113337300253E-002, & 4.802274679360077E-002, & 4.950592468304666E-002, & 5.083608261779847E-002, & 5.200910915174086E-002, & 5.302137852401040E-002, & 5.386976186571403E-002, & 5.455163687088878E-002, & 5.506489590176249E-002, & 5.540795250324512E-002, & 5.557974630651431E-002 /) ! Declare local variables character(len=80) :: imsl_file, in_file, out_file, var_name integer :: datid, fileid, i, ia, im1, ip1, j, ja, jm1, jp1, latdid, & latvid, londid, lonvid, mondid, month, monvid, status integer, dimension(nx, ny) :: imsl logical, dimension(nx, ny) :: masko real :: aoceana, aoceano, corr, error, sum, sumdy, total1, total2, weight real(kind=4), dimension(imt-2) :: lonts real(kind=4), dimension(jmt-2) :: latts real, dimension(ny) :: areaa real, dimension(imt) :: dxt real, dimension(jmt-2) :: areao real, dimension(jmt) :: cst, dyt, phi, phit real(kind=4), dimension(nx, ny, 12) :: aflux real(kind=4), dimension(imt-2, jmt-2, 12) :: ocean ! --------------------- ! BEGIN EXECUTABLE CODE ! --------------------- ! ---------------------------------------------------------------------- ! The following section calculates the values of CST, the cosine of the ! latitude of each tracer 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 cosine of the latitude of each tracer gridbox phi(1) = swldeg / radian phit(1) = phi(1) - 0.5 * dyt(1) / radius sumdy = phi(1) do j = 1, jmt if (j /= jmt) then sumdy = sumdy + dyt(j+1) / radius phi(j+1) = sumdy end if if (j /= 1) phit(j) = 0.5 * (phi(j-1) + phi(j)) cst(j) = cos(phit(j)) end do ! --------------------------------------------------------------------- ! Calculate DXT, the x-dimension of each gridbox, exactly as in ocean.f ! --------------------------------------------------------------------- ! Set the x-dimension of each box in degrees, and convert to centimetres do i = 2, imt-1 dxt(i) = 360.0 / float(imt-2) dxt(i) = dxt(i) * radius / radian end do ! Set cyclic conditions on DXT dxt(1) = dxt(imt-1) dxt(imt) = dxt(2) ! ------------------ ! Get the array IMSL ! ------------------ ! Get the name of the file containing IMSL call getarg(1, imsl_file) ! Get IMSL status = nf_open(imsl_file, nf_nowrite, fileid) status = nf_inq_varid(fileid, "imsl", datid) status = nf_get_var_int(fileid, datid, imsl) status = nf_close(fileid) ! ---------------------------------------------------- ! Initialise various arrays, exactly as in downscale.f ! ---------------------------------------------------- ! Calculate the areas of the AGCM gridboxes do j = 1, ny/2 areaa(j) = 2.0 * pi * w(j) * eradsq / float(nx) areaa(ny+1-j) = 2.0 * pi * w(j) * eradsq / float(nx) end do ! Calculate the areas of the OGCM gridboxes do j = 1, jmt-2 areao(j) = 1.0e-4 * cst(j+1) * dxt(1) * dyt(j+1) end do ! Create a land/sea mask do j = 1, ny do i = 1, nx if (imsl(i, j) .eq. 4) then masko(i, j) = .false. else masko(i, j) = .true. end if end do end do ! Calculate the surface area of the ocean according to the AGCM aoceana = 0.0 do j = 1, ny do i = 1, nx if (masko(i, j)) aoceana = aoceana + areaa(j) end do end do ! Calculate the surface area of the ocean according to the OGCM aoceano = 0.0 do j = 1, jmt-2 ja = (j + 1) / 2 do i = 1, imt-2 ia = (i + 1) / 2 if (masko(ia, ja)) aoceano = aoceano + areao(j) end do end do ! Display a summary write (*, *) write (*, *) "Ocean area according to AGCM = ", aoceana/1.0e6, " km^2" write (*, *) "Ocean area according to OGCM = ", aoceano/1.0e6, " km^2" write (*, *) ! --------------------------------------------------- ! Get the climatological surface fluxes from the AGCM ! --------------------------------------------------- ! Get the name of the input file call getarg(2, in_file) ! Get the variable name call getarg(3, var_name) ! Get the climatological surface fluxes from the AGCM status = nf_open(in_file, nf_nowrite, fileid) status = nf_inq_varid(fileid, var_name, datid) status = nf_get_var_real(fileid, datid, aflux) status = nf_close(fileid) ! --------------------------------------------------------------------------- ! The following section of code uses bilinear interpolation to downscale the ! surface flux fields from the AGCM grid to the new, high-resolution OGCM ! grid. ! ! These calculations are performed in exactly the same way as in downscale.f. ! --------------------------------------------------------------------------- ! Loop over months do month = 1, 12 ! Calculate the global integral of the AGCM surface fluxes total1 = 0.0 do j = 1, ny do i = 1, nx if (masko(i, j)) total1 = total1 + areaa(j) * aflux(i, j, month) end do end do ! Loop over all AGCM gridboxes do j = 1, ny jm1 = j - 1 jp1 = j + 1 do i = 1, nx im1 = i - 1 ip1 = i + 1 if (i .eq. 1) im1 = nx if (i .eq. nx) ip1 = 1 if (masko(i, j)) then ! This is an ocean gridbox, so use bilinear interpolation to ! calculate the surface flux on the OGCM grid ! ! (1) The OGCM gridbox at (2*i-1, 2*j-1) sum = 9.0 * areaa(j) * aflux(i, j, month) weight = 9.0 * areaa(j) if (masko(im1, j)) then sum = sum + 3.0 * areaa(j) * aflux(im1, j, month) weight = weight + 3.0 * areaa(j) end if if (j .gt. 1) then if (masko(i, jm1)) then sum = sum + 3.0 * areaa(jm1) * aflux(i, jm1, month) weight = weight + 3.0 * areaa(jm1) end if if (masko(im1, jm1)) then sum = sum + areaa(jm1) * aflux(im1, jm1, month) weight = weight + areaa(jm1) end if end if ocean(2*i-1, 2*j-1, month) = sum / weight ! (2) The OGCM gridbox at (2*i, 2*j-1) sum = 9.0 * areaa(j) * aflux(i, j, month) weight = 9.0 * areaa(j) if (masko(ip1, j)) then sum = sum + 3.0 * areaa(j) * aflux(ip1, j, month) weight = weight + 3.0 * areaa(j) end if if (j .gt. 1) then if (masko(i, jm1)) then sum = sum + 3.0 * areaa(jm1) * aflux(i, jm1, month) weight = weight + 3.0 * areaa(jm1) end if if (masko(ip1, jm1)) then sum = sum + areaa(jm1) * aflux(ip1, jm1, month) weight = weight + areaa(jm1) end if end if ocean(2*i, 2*j-1, month) = sum / weight ! (3) The OGCM gridbox at (2*i-1, 2*j) sum = 9.0 * areaa(j) * aflux(i, j, month) weight = 9.0 * areaa(j) if (masko(im1, j)) then sum = sum + 3.0 * areaa(j) * aflux(im1, j, month) weight = weight + 3.0 * areaa(j) end if if (j .lt. ny) then if (masko(i, jp1)) then sum = sum + 3.0 * areaa(jp1) * aflux(i, jp1, month) weight = weight + 3.0 * areaa(jp1) end if if (masko(im1, jp1)) then sum = sum + areaa(jp1) * aflux(im1, jp1, month) weight = weight + areaa(jp1) end if end if ocean(2*i-1, 2*j, month) = sum / weight ! (4) The OGCM gridbox at (2*i, 2*j) sum = 9.0 * areaa(j) * aflux(i, j, month) weight = 9.0 * areaa(j) if (masko(ip1, j)) then sum = sum + 3.0 * areaa(j) * aflux(ip1, j, month) weight = weight + 3.0 * areaa(j) end if if (j .lt. ny) then if (masko(i, jp1)) then sum = sum + 3.0 * areaa(jp1) * aflux(i, jp1, month) weight = weight + 3.0 * areaa(jp1) end if if (masko(ip1, jp1)) then sum = sum + areaa(jp1) * aflux(ip1, jp1, month) weight = weight + areaa(jp1) end if end if ocean(2*i, 2*j, month) = sum / weight else ! This is a land gridbox, so set the surface flux on the OGCM grid to ! zero ocean(2*i-1, 2*j-1, month) = 0.0 ocean(2*i, 2*j-1, month) = 0.0 ocean(2*i-1, 2*j, month) = 0.0 ocean(2*i, 2*j, month) = 0.0 end if end do end do ! Calculate the global integral of the OGCM surface fluxes total2 = 0.0 do j = 1, jmt-2 ja = (j + 1) / 2 do i = 1, imt-2 ia = (i + 1) / 2 if (masko(ia, ja)) total2 = total2 + areao(j) * ocean(i, j, month) end do end do ! Calculate the conservation error, and correct the OGCM surface fluxes error = total2 - total1 corr = error / aoceano do j = 1, jmt-2 ja = (j + 1) / 2 do i = 1, imt-2 ia = (i + 1) / 2 if (masko(ia, ja)) ocean(i, j, month) = ocean(i, j, month) - corr end do end do ! Display a summary write (*, *) write (*, *) "Month = ", month write (*, *) write (*, *) "TOTAL1 = ", total1 write (*, *) "TOTAL2 = ", total2 write (*, *) write (*, *) "ERROR = ", error write (*, *) write (*, *) "CORR = ", 0.0-corr write (*, *) end do ! -------------------------------------------------------- ! Derive the axis data for the interpolated surface fluxes ! -------------------------------------------------------- do i = 1, imt-2 lonts(i) = (real(i) - 1.5) * 360.0 / float(imt-2) end do do j = 1, jmt-2 latts(j) = phit(j+1) * radian end do ! -------------------------------- ! Save the data to the output file ! -------------------------------- ! Get the name of the output file call getarg(4, out_file) ! (1) Create netCDF file status = nf_create(out_file, nf_noclobber, fileid) ! (2) Define dimensions status = nf_def_dim(fileid, "lonts", imt-2, londid) status = nf_def_dim(fileid, "latts", jmt-2, latdid) status = nf_def_dim(fileid, "month", 12, mondid) ! (3) Define variables status = nf_def_var(fileid, "lonts", nf_float, 1, londid, lonvid) status = nf_def_var(fileid, "latts", nf_float, 1, latdid, latvid) status = nf_def_var(fileid, "month", nf_int, 1, mondid, monvid) status = nf_def_var(fileid, var_name, nf_float, 3, & (/ londid, latdid, mondid /), datid) ! (4) Put attributes 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, 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, lonts) status = nf_put_var_real(fileid, latvid, latts) 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, datid, ocean) ! (7) Close output file status = nf_close(fileid) end program downscale_athf