phylmd/Mobidic/regr_pr_av.f

module regr_pr_av_m

  implicit none

contains

  subroutine regr_pr_av(ncid, name, julien, paprs, v3)

    ! "regr_pr_av" stands for "regrid pressure averaging".

    ! This procedure reads a 2D latitude-pressure field from a NetCDF
    ! file, at a given day, regrids this field in pressure to the LMDZ
    ! vertical grid and packs it to the LMDZ horizontal "physics"
    ! grid. 

    ! We assume that, in the input file, the field has 3 dimensions:
    ! latitude, pressure, julian day.

    ! We assume that the input field is already on the LMDZ "rlatu"
    ! latitudes, except that latitudes are in ascending order in the
    ! input file.

    ! The target vertical LMDZ grid is the grid of layer boundaries.
    ! Regridding in pressure is done by averaging a step function of pressure.

    use dimens_m, only: iim, jjm, llm
    use dimphy, only: klon
    use grid_change, only: dyn_phy
    use netcdf95, only: nf95_inq_varid, nf95_get_var
    use nr_util, only: assert
    use numer_rec_95, only: regr1_step_av
    use press_coefoz_m, only: press_in_edg

    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):: 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
    ! layer "k".

    ! Variables local to the procedure:

    integer varid ! for NetCDF

    real v1(jjm + 1, size(press_in_edg) - 1)
    ! input field at day "julien"
    ! "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 on the "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, k

    !--------------------------------------------

    call assert(shape(v3) == (/klon, llm/), "regr_pr_av klon llm")
    call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs")

    call nf95_inq_varid(ncid, name, varid)

    ! Get data at the right day from the input file:
    call nf95_get_var(ncid, varid, v1, start=(/1, 1, julien/))
    ! Latitudes are in ascending order in the input file while
    ! "rlatu" is in descending order so we need to invert order:
    v1 = v1(jjm+1:1:-1, :)

    forall (k = 1:size(press_in_edg) - 1) v2(:, k) = pack(spread(v1(:, k), &
         dim = 1, ncopies = iim + 1), dyn_phy)

    ! 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

end module regr_pr_av_m
1	guez	90	module regr_pr_av_m
2	guez	3
3			implicit none
4
5			contains
6
7	guez	162	subroutine regr_pr_av(ncid, name, julien, paprs, v3)
8	guez	3
9	guez	90	! "regr_pr_av" stands for "regrid pressure averaging".
10
11			! This procedure reads a 2D latitude-pressure field from a NetCDF
12			! file, at a given day, regrids this field in pressure to the LMDZ
13			! vertical grid and packs it to the LMDZ horizontal "physics"
14			! grid.
15
16			! We assume that, in the input file, the field has 3 dimensions:
17			! latitude, pressure, julian day.
18
19			! We assume that the input field is already on the LMDZ "rlatu"
20			! latitudes, except that latitudes are in ascending order in the
21			! input file.
22
23	guez	19	! The target vertical LMDZ grid is the grid of layer boundaries.
24	guez	20	! Regridding in pressure is done by averaging a step function of pressure.
25	guez	3
26			use dimens_m, only: iim, jjm, llm
27			use dimphy, only: klon
28	guez	52	use grid_change, only: dyn_phy
29	guez	90	use netcdf95, only: nf95_inq_varid, nf95_get_var
30	guez	36	use nr_util, only: assert
31	guez	61	use numer_rec_95, only: regr1_step_av
32	guez	17	use press_coefoz_m, only: press_in_edg
33	guez	3
34			integer, intent(in):: ncid ! NetCDF ID of the file
35			character(len=*), intent(in):: name ! of the NetCDF variable
36	guez	7	integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
37	guez	3
38	guez	162	real, intent(in):: paprs(:, :) ! (klon, llm + 1)
39			! (pression pour chaque inter-couche, en Pa)
40
41	guez	7	real, intent(out):: v3(:, :) ! (klon, llm)
42	guez	90	! regridded field on the partial "physics" grid
43			! "v3(i, k)" is at longitude "xlon(i)", latitude "xlat(i)", in
44			! layer "k".
45	guez	7
46			! Variables local to the procedure:
47
48	guez	105	integer varid ! for NetCDF
49	guez	20
50	guez	162	real v1(jjm + 1, size(press_in_edg) - 1)
51	guez	90	! input field at day "julien"
52			! "v1(j, k)" is at latitude "rlatu(j)" and for
53			! pressure interval "[press_in_edg(k), press_in_edg(k+1)]".
54	guez	7
55	guez	162	real v2(klon, size(press_in_edg) - 1)
56			! Field on the "physics" horizontal grid. "v2(i, k)" is at
57			! longitude "xlon(i)", latitude "xlat(i)" and for pressure
58			! interval "[press_in_edg(k), press_in_edg(k+1)]".)
59	guez	7
60	guez	162	integer i, k
61	guez	20
62	guez	7	!--------------------------------------------
63
64	guez	162	call assert(shape(v3) == (/klon, llm/), "regr_pr_av klon llm")
65			call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs")
66	guez	7
67			call nf95_inq_varid(ncid, name, varid)
68
69			! Get data at the right day from the input file:
70	guez	90	call nf95_get_var(ncid, varid, v1, start=(/1, 1, julien/))
71	guez	20	! Latitudes are in ascending order in the input file while
72			! "rlatu" is in descending order so we need to invert order:
73	guez	7	v1 = v1(jjm+1:1:-1, :)
74
75	guez	162	forall (k = 1:size(press_in_edg) - 1) v2(:, k) = pack(spread(v1(:, k), &
76			dim = 1, ncopies = iim + 1), dyn_phy)
77
78	guez	7	! Regrid in pressure at each horizontal position:
79	guez	162	do i = 1, klon
80			v3(i, llm:1:-1) = regr1_step_av(v2(i, :), press_in_edg, &
81			paprs(i, llm+1:1:-1))
82			! (invert order of indices because "paprs" is in descending order)
83	guez	19	end do
84	guez	7
85	guez	90	end subroutine regr_pr_av
86	guez	7
87	guez	90	end module regr_pr_av_m