Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90 (revision 1158) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90 (revision 1158) @@ -0,0 +1,321 @@ +! $Id$ +module regr_lat_time_climoz_m + + ! Author: Lionel GUEZ + + implicit none + + private + public regr_lat_time_climoz + +contains + + subroutine regr_lat_time_climoz + + ! "regr_lat_time_climoz" stands for "regrid latitude time + ! climatology ozone". + + ! This procedure reads a climatology of ozone from a NetCDF file, + ! regrids it in latitude and time, and writes the regridded field + ! to a new NetCDF file. + + ! The input field depends on time, pressure level and latitude. + ! We assume that the input field is a step function of latitude + ! and that the input latitude coordinate gives the centers of steps. + ! Regridding in latitude is made by averaging, with a cosine of + ! latitude factor. + ! The target LMDZ latitude grid is the "scalar" grid: "rlatu". + ! The values of "rlatu" are taken to be the centers of intervals. + + ! Regridding in time is by linear interpolation. + ! Monthly values are processed to get daily values, on the basis + ! of a 360-day calendar. + ! The input file may contain either values for 12 months or values + ! for 14 months. + ! If there are 14 months then we assume that we have (in that order): + ! December, January, February, ..., November, December, January + ! We use the first December value to interpolate values between January + ! 1st and mid-January. + ! We use the last January value to interpolate values between + ! mid-December and end of December. + ! If there are only 12 months in the input file than we assume + ! periodicity for interpolation at the beginning and at the end of the + ! year. + + ! We assume that in the input file: + ! -- the latitude is in degrees and strictly monotonic (as all + ! NetCDF coordinate variables should be); + ! -- time increases (even though we do not use values of the input + ! time coordinate); + ! -- pressure is in hPa and in strictly ascending order (even + ! though we do not use pressure values here, we write the unit of + ! pressure in the NetCDF header, and we will use the assumptions later, + ! when we regrid in pressure). + + use regr1_step_av_m, only: regr1_step_av + use regr3_lint_m, only: regr3_lint + use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, & + nf95_put_var, nf95_gw_var + use netcdf, only: nf90_nowrite, nf90_get_var + + ! Variables local to the procedure: + + include "dimensions.h" + ! (for "jjm") + include "paramet.h" + ! (for the other included files) + include "comgeom2.h" + ! (for "rlatv") + include "comconst.h" + ! (for "pi") + + integer ncid_in, ncid_out ! NetCDF IDs for input and output files + integer n_plev ! number of pressure levels in the input data + integer n_lat! number of latitudes in the input data + + real, pointer:: latitude(:) + ! (of input data, converted to rad, sorted in strictly ascending order) + + real, allocatable:: lat_in_edg(:) + ! (edges of latitude intervals for input data, in rad, in strictly + ! ascending order) + + real, pointer:: plev(:) ! pressure level of input data + logical desc_lat ! latitude in descending order in the input file + + real, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13) + ! (ozone climatology from the input file) + ! ("o3_in(j, l, month)" is at latitude "latitude(j)" and pressure + ! level "plev(l)". "month" is between 1 and 12 or between 0 and 13) + + real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13) + ! (mean of "o3_in" over a latitude interval of LMDZ) + ! (First dimension is latitude interval. + ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)". + ! If "j" is between 2 and "jjm" then the interval is: + ! [rlatv(j), rlatv(j-1)] + ! If "j" is 1 or "jjm + 1" then the interval is: + ! [rlatv(1), pi / 2] + ! or: + ! [- pi / 2, rlatv(jjm)] + ! respectively. + ! "o3_regr_lat(:, l, :)" is for pressure level "plev(l)". + ! Last dimension is month number.) + + real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360) + ! (regridded ozone climatology) + ! ("o3_out(j, l, day)" is at latitude "rlatu(j)", pressure + ! level "plev(l)" and date "January 1st 0h" + "tmidday(day)", in a + ! 360-day calendar.) + + integer j + + integer varid_in, varid_out, varid_plev, varid_time + integer varid + ! (for NetCDF) + + real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * j, j = 0, 13)/) + ! (time to middle of month, in days since January 1st 0h, in a + ! 360-day calendar) + ! (We add values -15 and 375 so that, for example, day 3 of the year is + ! interpolated between the December and the January value.) + + real, parameter:: tmidday(360) = (/(j + 0.5, j = 0, 359)/) + ! (time to middle of day, in days since January 1st 0h, in a + ! 360-day calendar) + + !--------------------------------- + + print *, "Call sequence information: regr_lat_time_climoz" + + call nf95_open("climoz.nc", nf90_nowrite, ncid_in) + + ! Get coordinates from the input file: + + call nf95_inq_varid(ncid_in, "latitude", varid) + call nf95_gw_var(ncid_in, varid, latitude) + ! Convert from degrees to rad, because we will take the sine of latitude: + latitude = latitude / 180. * pi + n_lat = size(latitude) + ! We need to supply the latitudes to "regr1_step_av" in + ! ascending order, so invert order if necessary: + desc_lat = latitude(1) > latitude(n_lat) + if (desc_lat) latitude = latitude(n_lat:1:-1) + + ! Compute edges of latitude intervals: + allocate(lat_in_edg(n_lat + 1)) + lat_in_edg(1) = - pi / 2 + forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2 + lat_in_edg(n_lat + 1) = pi / 2 + deallocate(latitude) ! pointer + + call nf95_inq_varid(ncid_in, "plev", varid) + call nf95_gw_var(ncid_in, varid, plev) + n_plev = size(plev) + ! (We only need the pressure coordinate to copy it to the output file.) + + ! Create the output file and get the variable IDs: + call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & + varid_time) + + ! Write remaining coordinate variables: + call nf95_put_var(ncid_out, varid_plev, plev) + call nf95_put_var(ncid_out, varid_time, tmidday) + + deallocate(plev) ! pointer + + call nf95_inq_varid(ncid_in, "tro3", varid_in) + call nf95_gw_var(ncid_in, varid_in, o3_in) + if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :) + + call nf95_close(ncid_in) + + allocate(o3_regr_lat(jjm + 1, n_plev, 0:13)) + allocate(o3_out(jjm + 1, n_plev, 360)) + + ! Regrid in latitude: + ! We average with respect to sine of latitude, which is + ! equivalent to weighting by cosine of latitude: + if (size(o3_in, 3) == 12) then + print *, "Found 12 months in ozone climatology, assuming periodicity..." + o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, & + xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) + ! (invert order of indices in "o3_regr_lat" because "rlatu" is + ! in descending order) + + ! Duplicate January and December values, in preparation of time + ! interpolation: + o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12) + o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1) + else + print *, "Using 14 months in ozone climatology..." + o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, & + xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) + ! (invert order of indices in "o3_regr_lat" because "rlatu" is + ! in descending order) + end if + + deallocate(o3_in) ! pointer + + ! Regrid in time by linear interpolation: + o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday) + + ! Write to file: + call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :)) + ! (The order of "rlatu" is inverted in the output file) + + call nf95_close(ncid_out) + + end subroutine regr_lat_time_climoz + + !******************************************** + + subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & + varid_time) + + ! This subroutine creates the NetCDF output file, defines + ! dimensions and variables, and writes one of the coordinate variables. + + use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, & + nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var + use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global + + integer, intent(in):: ncid_in, n_plev + integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time + + ! Variables local to the procedure: + + include "dimensions.h" + ! (for "jjm") + include "paramet.h" + ! (for the other included files) + include "comgeom2.h" + ! (for "rlatu") + include "comconst.h" + ! (for "pi") + + integer ncerr + integer dimid_rlatu, dimid_plev, dimid_time + integer varid_rlatu + + !--------------------------- + + print *, "Call sequence information: prepare_out" + + call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out) + + ! Dimensions: + call nf95_def_dim(ncid_out, "time", 360, dimid_time) + call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev) + call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu) + + ! Define coordinate variables: + + call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time) + call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1") + call nf95_put_att(ncid_out, varid_time, "calendar", "360_day") + call nf95_put_att(ncid_out, varid_time, "standard_name", "time") + + call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev) + call nf95_put_att(ncid_out, varid_plev, "units", "millibar") + call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure") + call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure") + + call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu) + call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north") + call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude") + + ! Define the primary variable: + + call nf95_def_var(ncid_out, "tro3", nf90_float, & + (/dimid_rlatu, dimid_plev, dimid_time/), varid_out) + call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction") + call nf95_put_att(ncid_out, varid_out, "standard_name", & + "mole_fraction_of_ozone_in_air") + + ! Global attributes: + + call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ") + + ! The following commands, copying attributes, may fail. + ! That is OK. + ! It should just mean that the attribute is not defined in the input file. + + call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, & + nf90_global, ncerr) + call handle_err_copy_att("Conventions") + + call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, & + ncerr) + call handle_err_copy_att("title") + + call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, & + ncerr) + call handle_err_copy_att("source") + + call nf95_enddef(ncid_out) + + ! Write one of the coordinate variables: + call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.) + ! (convert from rad to degrees and sort in ascending order) + + contains + + subroutine handle_err_copy_att(att_name) + + use netcdf, only: nf90_noerr, nf90_strerror + + character(len=*), intent(in):: att_name + + !---------------------------------------- + + if (ncerr /= nf90_noerr) then + print *, "regr_lat_time_climoz_m prepare_out nf90_copy_att " & + // att_name // " -- " // trim(nf90_strerror(ncerr)) + end if + + end subroutine handle_err_copy_att + + end subroutine prepare_out + +end module regr_lat_time_climoz_m Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.f90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.f90 (revision 1157) +++ (revision ) @@ -1,321 +1,0 @@ -! $Id$ -module regr_lat_time_climoz_m - - ! Author: Lionel GUEZ - - implicit none - - private - public regr_lat_time_climoz - -contains - - subroutine regr_lat_time_climoz - - ! "regr_lat_time_climoz" stands for "regrid latitude time - ! climatology ozone". - - ! This procedure reads a climatology of ozone from a NetCDF file, - ! regrids it in latitude and time, and writes the regridded field - ! to a new NetCDF file. - - ! The input field depends on time, pressure level and latitude. - ! We assume that the input field is a step function of latitude - ! and that the input latitude coordinate gives the centers of steps. - ! Regridding in latitude is made by averaging, with a cosine of - ! latitude factor. - ! The target LMDZ latitude grid is the "scalar" grid: "rlatu". - ! The values of "rlatu" are taken to be the centers of intervals. - - ! Regridding in time is by linear interpolation. - ! Monthly values are processed to get daily values, on the basis - ! of a 360-day calendar. - ! The input file may contain either values for 12 months or values - ! for 14 months. - ! If there are 14 months then we assume that we have (in that order): - ! December, January, February, ..., November, December, January - ! We use the first December value to interpolate values between January - ! 1st and mid-January. - ! We use the last January value to interpolate values between - ! mid-December and end of December. - ! If there are only 12 months in the input file than we assume - ! periodicity for interpolation at the beginning and at the end of the - ! year. - - ! We assume that in the input file: - ! -- the latitude is in degrees and strictly monotonic (as all - ! NetCDF coordinate variables should be); - ! -- time increases (even though we do not use values of the input - ! time coordinate); - ! -- pressure is in hPa and in strictly ascending order (even - ! though we do not use pressure values here, we write the unit of - ! pressure in the NetCDF header, and we will use the assumptions later, - ! when we regrid in pressure). - - use regr1_step_av_m, only: regr1_step_av - use regr3_lint_m, only: regr3_lint - use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, & - nf95_put_var, nf95_gw_var - use netcdf, only: nf90_nowrite, nf90_get_var - - ! Variables local to the procedure: - - include "dimensions.h" - ! (for "jjm") - include "paramet.h" - ! (for the other included files) - include "comgeom2.h" - ! (for "rlatv") - include "comconst.h" - ! (for "pi") - - integer ncid_in, ncid_out ! NetCDF IDs for input and output files - integer n_plev ! number of pressure levels in the input data - integer n_lat! number of latitudes in the input data - - real, pointer:: latitude(:) - ! (of input data, converted to rad, sorted in strictly ascending order) - - real, allocatable:: lat_in_edg(:) - ! (edges of latitude intervals for input data, in rad, in strictly - ! ascending order) - - real, pointer:: plev(:) ! pressure level of input data - logical desc_lat ! latitude in descending order in the input file - - real, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13) - ! (ozone climatology from the input file) - ! ("o3_in(j, l, month)" is at latitude "latitude(j)" and pressure - ! level "plev(l)". "month" is between 1 and 12 or between 0 and 13) - - real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13) - ! (mean of "o3_in" over a latitude interval of LMDZ) - ! (First dimension is latitude interval. - ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)". - ! If "j" is between 2 and "jjm" then the interval is: - ! [rlatv(j), rlatv(j-1)] - ! If "j" is 1 or "jjm + 1" then the interval is: - ! [rlatv(1), pi / 2] - ! or: - ! [- pi / 2, rlatv(jjm)] - ! respectively. - ! "o3_regr_lat(:, l, :)" is for pressure level "plev(l)". - ! Last dimension is month number.) - - real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360) - ! (regridded ozone climatology) - ! ("o3_out(j, l, day)" is at latitude "rlatu(j)", pressure - ! level "plev(l)" and date "January 1st 0h" + "tmidday(day)", in a - ! 360-day calendar.) - - integer j - - integer varid_in, varid_out, varid_plev, varid_time - integer varid - ! (for NetCDF) - - real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * j, j = 0, 13)/) - ! (time to middle of month, in days since January 1st 0h, in a - ! 360-day calendar) - ! (We add values -15 and 375 so that, for example, day 3 of the year is - ! interpolated between the December and the January value.) - - real, parameter:: tmidday(360) = (/(j + 0.5, j = 0, 359)/) - ! (time to middle of day, in days since January 1st 0h, in a - ! 360-day calendar) - - !--------------------------------- - - print *, "Call sequence information: regr_lat_time_climoz" - - call nf95_open("climoz.nc", nf90_nowrite, ncid_in) - - ! Get coordinates from the input file: - - call nf95_inq_varid(ncid_in, "latitude", varid) - call nf95_gw_var(ncid_in, varid, latitude) - ! Convert from degrees to rad, because we will take the sine of latitude: - latitude = latitude / 180. * pi - n_lat = size(latitude) - ! We need to supply the latitudes to "regr1_step_av" in - ! ascending order, so invert order if necessary: - desc_lat = latitude(1) > latitude(n_lat) - if (desc_lat) latitude = latitude(n_lat:1:-1) - - ! Compute edges of latitude intervals: - allocate(lat_in_edg(n_lat + 1)) - lat_in_edg(1) = - pi / 2 - forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2 - lat_in_edg(n_lat + 1) = pi / 2 - deallocate(latitude) ! pointer - - call nf95_inq_varid(ncid_in, "plev", varid) - call nf95_gw_var(ncid_in, varid, plev) - n_plev = size(plev) - ! (We only need the pressure coordinate to copy it to the output file.) - - ! Create the output file and get the variable IDs: - call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & - varid_time) - - ! Write remaining coordinate variables: - call nf95_put_var(ncid_out, varid_plev, plev) - call nf95_put_var(ncid_out, varid_time, tmidday) - - deallocate(plev) ! pointer - - call nf95_inq_varid(ncid_in, "tro3", varid_in) - call nf95_gw_var(ncid_in, varid_in, o3_in) - if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :) - - call nf95_close(ncid_in) - - allocate(o3_regr_lat(jjm + 1, n_plev, 0:13)) - allocate(o3_out(jjm + 1, n_plev, 360)) - - ! Regrid in latitude: - ! We average with respect to sine of latitude, which is - ! equivalent to weighting by cosine of latitude: - if (size(o3_in, 3) == 12) then - print *, "Found 12 months in ozone climatology, assuming periodicity..." - o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, & - xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) - ! (invert order of indices in "o3_regr_lat" because "rlatu" is - ! in descending order) - - ! Duplicate January and December values, in preparation of time - ! interpolation: - o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12) - o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1) - else - print *, "Using 14 months in ozone climatology..." - o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, & - xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) - ! (invert order of indices in "o3_regr_lat" because "rlatu" is - ! in descending order) - end if - - deallocate(o3_in) ! pointer - - ! Regrid in time by linear interpolation: - o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday) - - ! Write to file: - call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :)) - ! (The order of "rlatu" is inverted in the output file) - - call nf95_close(ncid_out) - - end subroutine regr_lat_time_climoz - - !******************************************** - - subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & - varid_time) - - ! This subroutine creates the NetCDF output file, defines - ! dimensions and variables, and writes one of the coordinate variables. - - use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, & - nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var - use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global - - integer, intent(in):: ncid_in, n_plev - integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time - - ! Variables local to the procedure: - - include "dimensions.h" - ! (for "jjm") - include "paramet.h" - ! (for the other included files) - include "comgeom2.h" - ! (for "rlatu") - include "comconst.h" - ! (for "pi") - - integer ncerr - integer dimid_rlatu, dimid_plev, dimid_time - integer varid_rlatu - - !--------------------------- - - print *, "Call sequence information: prepare_out" - - call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out) - - ! Dimensions: - call nf95_def_dim(ncid_out, "time", 360, dimid_time) - call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev) - call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu) - - ! Define coordinate variables: - - call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time) - call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1") - call nf95_put_att(ncid_out, varid_time, "calendar", "360_day") - call nf95_put_att(ncid_out, varid_time, "standard_name", "time") - - call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev) - call nf95_put_att(ncid_out, varid_plev, "units", "millibar") - call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure") - call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure") - - call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu) - call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north") - call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude") - - ! Define the primary variable: - - call nf95_def_var(ncid_out, "tro3", nf90_float, & - (/dimid_rlatu, dimid_plev, dimid_time/), varid_out) - call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction") - call nf95_put_att(ncid_out, varid_out, "standard_name", & - "mole_fraction_of_ozone_in_air") - - ! Global attributes: - - call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ") - - ! The following commands, copying attributes, may fail. - ! That is OK. - ! It should just mean that the attribute is not defined in the input file. - - call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, & - nf90_global, ncerr) - call handle_err_copy_att("Conventions") - - call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, & - ncerr) - call handle_err_copy_att("title") - - call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, & - ncerr) - call handle_err_copy_att("source") - - call nf95_enddef(ncid_out) - - ! Write one of the coordinate variables: - call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.) - ! (convert from rad to degrees and sort in ascending order) - - contains - - subroutine handle_err_copy_att(att_name) - - use netcdf, only: nf90_noerr, nf90_strerror - - character(len=*), intent(in):: att_name - - !---------------------------------------- - - if (ncerr /= nf90_noerr) then - print *, "regr_lat_time_climoz_m prepare_out nf90_copy_att " & - // att_name // " -- " // trim(nf90_strerror(ncerr)) - end if - - end subroutine handle_err_copy_att - - end subroutine prepare_out - -end module regr_lat_time_climoz_m Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.F90 (revision 1158) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.F90 (revision 1158) @@ -0,0 +1,191 @@ +! $Id$ +module regr_pr + + ! In both procedures of this module: + ! -- the root process reads a 2D latitude-pressure field from a + ! NetCDF file, at a given day. + ! -- the field is packed to the LMDZ horizontal "physics" + ! grid and scattered to all threads of all processes; + ! -- in all the threads of all the processes, the field is regridded in + ! pressure to the LMDZ vertical grid. + ! We assume that, in the input file, the field has 3 dimensions: + ! latitude, pressure, julian day. + ! We assume that latitudes are in ascending order in the input file. + + + implicit none + +contains + + subroutine regr_pr_av(ncid, name, julien, press_in_edg, paprs, v3) + + ! "regr_pr_av" stands for "regrid pressure averaging". + ! The target vertical LMDZ grid is the grid of layer boundaries. + ! Regridding in pressure is done by averaging a step function of pressure. + + use dimphy, only: klon + use netcdf95, only: nf95_inq_varid, handle_err + use netcdf, only: nf90_get_var + use assert_m, only: assert + use regr1_step_av_m, only: regr1_step_av + use mod_phys_lmdz_mpi_data, only: is_mpi_root + + use mod_phys_lmdz_transfert_para, only: scatter2d + ! (pack to the LMDZ horizontal "physics" grid and scatter) + + integer, intent(in):: ncid ! NetCDF ID of the file + character(len=*), intent(in):: name ! of the NetCDF variable + integer, intent(in):: julien ! jour julien, 1 <= julien <= 360 + + real, intent(in):: press_in_edg(:) + ! (edges of pressure intervals for input data, in Pa, in strictly + ! ascending order) + + real, intent(in):: paprs(:, :) ! (klon, llm + 1) + ! (pression pour chaque inter-couche, en Pa) + + real, intent(out):: v3(:, :) ! (klon, llm) + ! (regridded field on the partial "physics" grid) + ! ("v3(i, k)" is at longitude "xlon(i)", latitude + ! "xlat(i)", in pressure interval "[paprs(i, l+1), paprs(i, l)]".) + + ! Variables local to the procedure: + + include "dimensions.h" + integer varid, ncerr ! for NetCDF + + real v1(iim, jjm + 1, size(press_in_edg) - 1) + ! (input field at day "julien", on the global "dynamics" horizontal grid) + ! (First dimension is for longitude. + ! The value is the same for all longitudes. + ! "v1(:, j, k)" is at latitude "rlatu(j)" and for + ! pressure interval "[press_in_edg(k), press_in_edg(k+1)]".) + + real v2(klon, size(press_in_edg) - 1) + ! (field scattered to the partial "physics" horizontal grid) + ! ("v2(i, k)" is at longitude "xlon(i)", latitude "xlat(i)" and + ! for pressure interval "[press_in_edg(k), press_in_edg(k+1)]".) + + integer i + + !-------------------------------------------- + + call assert(shape(v3) == (/klon, llm/), "regr_pr_av v3") + call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs") + + !$omp master + if (is_mpi_root) then + call nf95_inq_varid(ncid, name, varid) + + ! Get data at the right day from the input file: + ncerr = nf90_get_var(ncid, varid, v1(1, :, :), start=(/1, 1, julien/)) + call handle_err("regr_pr_av nf90_get_var " // name, ncerr, ncid) + ! Latitudes are in ascending order in the input file while + ! "rlatu" is in descending order so we need to invert order: + v1(1, :, :) = v1(1, jjm+1:1:-1, :) + + ! Duplicate on all longitudes: + v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1) + end if + !$omp end master + + call scatter2d(v1, v2) + + ! Regrid in pressure at each horizontal position: + do i = 1, klon + v3(i, llm:1:-1) = regr1_step_av(v2(i, :), press_in_edg, & + paprs(i, llm+1:1:-1)) + ! (invert order of indices because "paprs" is in descending order) + end do + + end subroutine regr_pr_av + + !*************************************************************** + + subroutine regr_pr_int(ncid, name, julien, plev, pplay, top_value, v3) + + ! "regr_pr_int" stands for "regrid pressure interpolation". + ! The target vertical LMDZ grid is the grid of mid-layers. + ! Regridding is by linear interpolation. + + use dimphy, only: klon + use netcdf95, only: nf95_inq_varid, handle_err + use netcdf, only: nf90_get_var + use assert_m, only: assert + use regr1_lint_m, only: regr1_lint + use mod_phys_lmdz_mpi_data, only: is_mpi_root + + use mod_phys_lmdz_transfert_para, only: scatter2d + ! (pack to the LMDZ horizontal "physics" grid and scatter) + + integer, intent(in):: ncid ! NetCDF ID of the file + character(len=*), intent(in):: name ! of the NetCDF variable + integer, intent(in):: julien ! jour julien, 1 <= julien <= 360 + + real, intent(in):: plev(:) + ! (pressure level of input data, in Pa, in strictly ascending order) + + real, intent(in):: pplay(:, :) ! (klon, llm) + ! (pression pour le mileu de chaque couche, en Pa) + + real, intent(in):: top_value + ! (extra value of field at 0 pressure) + + real, intent(out):: v3(:, :) ! (klon, llm) + ! (regridded field on the partial "physics" grid) + ! ("v3(i, k)" is at longitude "xlon(i)", latitude + ! "xlat(i)", middle of layer "k".) + + ! Variables local to the procedure: + + include "dimensions.h" + integer varid, ncerr ! for NetCDF + + real v1(iim, jjm + 1, 0:size(plev)) + ! (input field at day "julien", on the global "dynamics" horizontal grid) + ! (First dimension is for longitude. + ! The value is the same for all longitudes. + ! "v1(:, j, k >=1)" is at latitude "rlatu(j)" and pressure "plev(k)".) + + real v2(klon, 0:size(plev)) + ! (field scattered to the partial "physics" horizontal grid) + ! "v2(i, k >= 1)" is at longitude "xlon(i)", latitude "xlat(i)" + ! and pressure "plev(k)".) + + integer i + + !-------------------------------------------- + + call assert(shape(v3) == (/klon, llm/), "regr_pr_int v3") + call assert(shape(pplay) == (/klon, llm/), "regr_pr_int pplay") + + !$omp master + if (is_mpi_root) then + call nf95_inq_varid(ncid, name, varid) + + ! Get data at the right day from the input file: + ncerr = nf90_get_var(ncid, varid, v1(1, :, 1:), start=(/1, 1, julien/)) + call handle_err("regr_pr_int nf90_get_var " // name, ncerr, ncid) + ! Latitudes are in ascending order in the input file while + ! "rlatu" is in descending order so we need to invert order: + v1(1, :, 1:) = v1(1, jjm+1:1:-1, 1:) + + ! Complete "v1" with the value at 0 pressure: + v1(1, :, 0) = top_value + + ! Duplicate on all longitudes: + v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1) + end if + !$omp end master + + call scatter2d(v1, v2) + + ! Regrid in pressure at each horizontal position: + do i = 1, klon + v3(i, llm:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, llm:1:-1)) + ! (invert order of indices because "pplay" is in descending order) + end do + + end subroutine regr_pr_int + +end module regr_pr Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.f90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.f90 (revision 1157) +++ (revision ) @@ -1,191 +1,0 @@ -! $Id$ -module regr_pr - - ! In both procedures of this module: - ! -- the root process reads a 2D latitude-pressure field from a - ! NetCDF file, at a given day. - ! -- the field is packed to the LMDZ horizontal "physics" - ! grid and scattered to all threads of all processes; - ! -- in all the threads of all the processes, the field is regridded in - ! pressure to the LMDZ vertical grid. - ! We assume that, in the input file, the field has 3 dimensions: - ! latitude, pressure, julian day. - ! We assume that latitudes are in ascending order in the input file. - - - implicit none - -contains - - subroutine regr_pr_av(ncid, name, julien, press_in_edg, paprs, v3) - - ! "regr_pr_av" stands for "regrid pressure averaging". - ! The target vertical LMDZ grid is the grid of layer boundaries. - ! Regridding in pressure is done by averaging a step function of pressure. - - use dimphy, only: klon - use netcdf95, only: nf95_inq_varid, handle_err - use netcdf, only: nf90_get_var - use assert_m, only: assert - use regr1_step_av_m, only: regr1_step_av - use mod_phys_lmdz_mpi_data, only: is_mpi_root - - use mod_phys_lmdz_transfert_para, only: scatter2d - ! (pack to the LMDZ horizontal "physics" grid and scatter) - - integer, intent(in):: ncid ! NetCDF ID of the file - character(len=*), intent(in):: name ! of the NetCDF variable - integer, intent(in):: julien ! jour julien, 1 <= julien <= 360 - - real, intent(in):: press_in_edg(:) - ! (edges of pressure intervals for input data, in Pa, in strictly - ! ascending order) - - real, intent(in):: paprs(:, :) ! (klon, llm + 1) - ! (pression pour chaque inter-couche, en Pa) - - real, intent(out):: v3(:, :) ! (klon, llm) - ! (regridded field on the partial "physics" grid) - ! ("v3(i, k)" is at longitude "xlon(i)", latitude - ! "xlat(i)", in pressure interval "[paprs(i, l+1), paprs(i, l)]".) - - ! Variables local to the procedure: - - include "dimensions.h" - integer varid, ncerr ! for NetCDF - - real v1(iim, jjm + 1, size(press_in_edg) - 1) - ! (input field at day "julien", on the global "dynamics" horizontal grid) - ! (First dimension is for longitude. - ! The value is the same for all longitudes. - ! "v1(:, j, k)" is at latitude "rlatu(j)" and for - ! pressure interval "[press_in_edg(k), press_in_edg(k+1)]".) - - real v2(klon, size(press_in_edg) - 1) - ! (field scattered to the partial "physics" horizontal grid) - ! ("v2(i, k)" is at longitude "xlon(i)", latitude "xlat(i)" and - ! for pressure interval "[press_in_edg(k), press_in_edg(k+1)]".) - - integer i - - !-------------------------------------------- - - call assert(shape(v3) == (/klon, llm/), "regr_pr_av v3") - call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs") - - !$omp master - if (is_mpi_root) then - call nf95_inq_varid(ncid, name, varid) - - ! Get data at the right day from the input file: - ncerr = nf90_get_var(ncid, varid, v1(1, :, :), start=(/1, 1, julien/)) - call handle_err("regr_pr_av nf90_get_var " // name, ncerr, ncid) - ! Latitudes are in ascending order in the input file while - ! "rlatu" is in descending order so we need to invert order: - v1(1, :, :) = v1(1, jjm+1:1:-1, :) - - ! Duplicate on all longitudes: - v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1) - end if - !$omp end master - - call scatter2d(v1, v2) - - ! Regrid in pressure at each horizontal position: - do i = 1, klon - v3(i, llm:1:-1) = regr1_step_av(v2(i, :), press_in_edg, & - paprs(i, llm+1:1:-1)) - ! (invert order of indices because "paprs" is in descending order) - end do - - end subroutine regr_pr_av - - !*************************************************************** - - subroutine regr_pr_int(ncid, name, julien, plev, pplay, top_value, v3) - - ! "regr_pr_int" stands for "regrid pressure interpolation". - ! The target vertical LMDZ grid is the grid of mid-layers. - ! Regridding is by linear interpolation. - - use dimphy, only: klon - use netcdf95, only: nf95_inq_varid, handle_err - use netcdf, only: nf90_get_var - use assert_m, only: assert - use regr1_lint_m, only: regr1_lint - use mod_phys_lmdz_mpi_data, only: is_mpi_root - - use mod_phys_lmdz_transfert_para, only: scatter2d - ! (pack to the LMDZ horizontal "physics" grid and scatter) - - integer, intent(in):: ncid ! NetCDF ID of the file - character(len=*), intent(in):: name ! of the NetCDF variable - integer, intent(in):: julien ! jour julien, 1 <= julien <= 360 - - real, intent(in):: plev(:) - ! (pressure level of input data, in Pa, in strictly ascending order) - - real, intent(in):: pplay(:, :) ! (klon, llm) - ! (pression pour le mileu de chaque couche, en Pa) - - real, intent(in):: top_value - ! (extra value of field at 0 pressure) - - real, intent(out):: v3(:, :) ! (klon, llm) - ! (regridded field on the partial "physics" grid) - ! ("v3(i, k)" is at longitude "xlon(i)", latitude - ! "xlat(i)", middle of layer "k".) - - ! Variables local to the procedure: - - include "dimensions.h" - integer varid, ncerr ! for NetCDF - - real v1(iim, jjm + 1, 0:size(plev)) - ! (input field at day "julien", on the global "dynamics" horizontal grid) - ! (First dimension is for longitude. - ! The value is the same for all longitudes. - ! "v1(:, j, k >=1)" is at latitude "rlatu(j)" and pressure "plev(k)".) - - real v2(klon, 0:size(plev)) - ! (field scattered to the partial "physics" horizontal grid) - ! "v2(i, k >= 1)" is at longitude "xlon(i)", latitude "xlat(i)" - ! and pressure "plev(k)".) - - integer i - - !-------------------------------------------- - - call assert(shape(v3) == (/klon, llm/), "regr_pr_int v3") - call assert(shape(pplay) == (/klon, llm/), "regr_pr_int pplay") - - !$omp master - if (is_mpi_root) then - call nf95_inq_varid(ncid, name, varid) - - ! Get data at the right day from the input file: - ncerr = nf90_get_var(ncid, varid, v1(1, :, 1:), start=(/1, 1, julien/)) - call handle_err("regr_pr_int nf90_get_var " // name, ncerr, ncid) - ! Latitudes are in ascending order in the input file while - ! "rlatu" is in descending order so we need to invert order: - v1(1, :, 1:) = v1(1, jjm+1:1:-1, 1:) - - ! Complete "v1" with the value at 0 pressure: - v1(1, :, 0) = top_value - - ! Duplicate on all longitudes: - v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1) - end if - !$omp end master - - call scatter2d(v1, v2) - - ! Regrid in pressure at each horizontal position: - do i = 1, klon - v3(i, llm:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, llm:1:-1)) - ! (invert order of indices because "pplay" is in descending order) - end do - - end subroutine regr_pr_int - -end module regr_pr