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	module regr_pr_av_m
2
3	implicit none
4
5	contains
6
7	subroutine regr_pr_av(ncid, name, julien, paprs, v3)
8
9	! "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	! The target vertical LMDZ grid is the grid of layer boundaries.
24	! Regridding in pressure is done by averaging a step function of pressure.
25
26	use dimens_m, only: iim, jjm, llm
27	use dimphy, only: klon
28	use grid_change, only: dyn_phy
29	use netcdf95, only: nf95_inq_varid, nf95_get_var
30	use nr_util, only: assert
31	use numer_rec_95, only: regr1_step_av
32	use press_coefoz_m, only: press_in_edg
33
34	integer, intent(in):: ncid ! NetCDF ID of the file
35	character(len=*), intent(in):: name ! of the NetCDF variable
36	integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
37
38	real, intent(in):: paprs(:, :) ! (klon, llm + 1)
39	! (pression pour chaque inter-couche, en Pa)
40
41	real, intent(out):: v3(:, :) ! (klon, llm)
42	! regridded field on the partial "physics" grid
43	! "v3(i, k)" is at longitude "xlon(i)", latitude "xlat(i)", in
44	! layer "k".
45
46	! Variables local to the procedure:
47
48	integer varid ! for NetCDF
49
50	real v1(jjm + 1, size(press_in_edg) - 1)
51	! 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
55	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
60	integer i, k
61
62	!--------------------------------------------
63
64	call assert(shape(v3) == (/klon, llm/), "regr_pr_av klon llm")
65	call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs")
66
67	call nf95_inq_varid(ncid, name, varid)
68
69	! Get data at the right day from the input file:
70	call nf95_get_var(ncid, varid, v1, start=(/1, 1, julien/))
71	! Latitudes are in ascending order in the input file while
72	! "rlatu" is in descending order so we need to invert order:
73	v1 = v1(jjm+1:1:-1, :)
74
75	forall (k = 1:size(press_in_edg) - 1) v2(:, k) = pack(spread(v1(:, k), &
76	dim = 1, ncopies = iim + 1), dyn_phy)
77
78	! Regrid in pressure at each horizontal position:
79	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	end do
84
85	end subroutine regr_pr_av
86
87	end module regr_pr_av_m