libf/dyn3d/startdyn.f90

MODULE startdyn

  ! From startvar.F, version 1.4
  ! 2006/01/27 15:14:22 Fairhead

  IMPLICIT NONE

  private
  public start_init_dyn, start_inter_3d

  INTEGER fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn

  REAL, ALLOCATABLE:: lon_ini(:), lat_ini(:)
  ! (longitude and latitude from the input file, converted to rad)

  real, ALLOCATABLE:: levdyn_ini(:)

CONTAINS

  SUBROUTINE start_init_dyn(tsol_2d, psol)

    ! Host associated variables appearing and modified in this procedure :
    ! iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn, lon_ini, lat_ini, levdyn_ini

    USE flincom, only: flininfo, flinopen_nozoom
    use comgeom, only: aire_2d, apoln, apols
    use conf_dat2d_m, only: conf_dat2d
    use inter_barxy_m, only: inter_barxy
    use comgeom, only: rlonu, rlatv
    use dimens_m, only: iim, jjm
    use gr_int_dyn_m, only: gr_int_dyn
    use start_init_orog_m, only: phis
    use start_init_phys_m, only: start_init_phys
    use nr_util, only: assert, pi
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid

    REAL, intent(in):: tsol_2d(:, :) ! (iim + 1, jjm + 1)
    REAL, intent(out):: psol(:, :) ! (iim + 1, jjm + 1) surface pressure, in Pa

    ! Local:

    REAL date, dt
    INTEGER itau(1), ncid, varid
    REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:)

    REAL, ALLOCATABLE:: lon_dyn(:, :), lat_dyn(:, :)
    ! (longitude and latitude from the input file, in rad or degrees)

    REAL, ALLOCATABLE:: var_ana(:, :)
    real z(iim + 1, jjm + 1)
    real tmp_var(iim, jjm + 1)

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

    print *, "Call sequence information: start_init_dyn"
    call assert((/size(tsol_2d, 1), size(psol, 1)/) == iim + 1, &
         "start_init_dyn size 1")
    call assert((/size(tsol_2d, 2), size(psol, 2)/) == jjm + 1, &
         "start_init_dyn size 2")
    CALL flininfo('ECDYN.nc', iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
    print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &
         ", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn

    ALLOCATE(lat_dyn(iml_dyn, jml_dyn))
    ALLOCATE(lon_dyn(iml_dyn, jml_dyn))
    ALLOCATE(levdyn_ini(llm_dyn))

    CALL flinopen_nozoom(iml_dyn, jml_dyn, llm_dyn, &
         lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)

    ALLOCATE(var_ana(iml_dyn, jml_dyn))
    ALLOCATE(lon_rad(iml_dyn))
    ALLOCATE(lon_ini(iml_dyn))

    IF (MAXVAL(lon_dyn) > pi) THEN
       ! Assume "lon_dyn" is in degrees 
       lon_ini = lon_dyn(:, 1) * pi / 180.
    ELSE
       lon_ini = lon_dyn(:, 1) 
    ENDIF

    ALLOCATE(lat_rad(jml_dyn))
    ALLOCATE(lat_ini(jml_dyn))

    IF (MAXVAL(lat_dyn) > pi) THEN
       lat_ini = lat_dyn(1, :) * pi / 180.
    ELSE
       lat_ini = lat_dyn(1, :) 
    ENDIF

    call nf95_open('ECDYN.nc', nf90_nowrite, ncid)

    ! 'Z': Surface geopotential
    call nf95_inq_varid(ncid, 'Z', varid)
    call nf95_get_var(ncid, varid, var_ana)
    CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
    CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
         rlatv, tmp_var)
    z = gr_int_dyn(tmp_var)

    ! 'SP': Surface pressure
    call nf95_inq_varid(ncid, 'SP', varid)
    call nf95_get_var(ncid, varid, var_ana)
    CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
    CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
         rlatv, tmp_var) 
    psol = gr_int_dyn(tmp_var)

    call nf95_close(ncid)

    psol(:iim, :) = psol(:iim, :) &
         * (1. + (z(:iim, :) - phis(:iim, :)) / 287. / tsol_2d(:iim, :))
    psol(iim + 1, :) = psol(1, :)

    psol(:, 1) = SUM(aire_2d(:iim, 1) * psol(:iim, 1)) / apoln
    psol(:, jjm + 1) = SUM(aire_2d(:iim, jjm + 1) * psol(:iim, jjm + 1)) &
         / apols

  END SUBROUTINE start_init_dyn

  !********************************

  subroutine start_inter_3d(varname, lon_in2, lat_in2, pls_in, var3d)

    ! This procedure gets a 3D variable from a file and interpolates it.

    use nr_util, only: assert_eq
    use numer_rec, only: spline, splint
    use inter_barxy_m, only: inter_barxy
    use gr_int_dyn_m, only: gr_int_dyn
    use conf_dat3d_m, only: conf_dat3d
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid

    CHARACTER(len=*), intent(in):: varname
    REAL, intent(in):: lon_in2(:) ! (iml)
    REAL, intent(in):: lat_in2(:)
    REAL, intent(in):: pls_in(:, :, :) ! (iml, jml, lml)
    REAL, intent(out):: var3d(:, :, :) ! (iml, jml, lml)

    ! LOCAL:
    INTEGER iml, jml, lml, ncid, varid
    INTEGER ii, ij, il
    REAL lon_rad(iml_dyn), lat_rad(jml_dyn)
    REAL lev_dyn(llm_dyn)
    REAL var_tmp2d(size(lon_in2)-1, size(pls_in, 2))
    real var_tmp3d(size(lon_in2), size(pls_in, 2), llm_dyn)
    REAL ax(llm_dyn), ay(llm_dyn), yder(llm_dyn)
    real var_ana3d(iml_dyn, jml_dyn, llm_dyn)

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

    print *, "Call sequence information: start_inter_3d"

    iml = assert_eq(size(pls_in, 1), size(lon_in2), size(var3d, 1), &
         "start_inter_3d iml")
    jml = assert_eq(size(pls_in, 2), size(var3d, 2), "start_inter_3d jml")
    lml = assert_eq(size(pls_in, 3), size(var3d, 3), "start_inter_3d lml")

    print *, "iml = ", iml, ", jml = ", jml
    print *, "varname = ", varname
    call nf95_open('ECDYN.nc', nf90_nowrite, ncid)
    call nf95_inq_varid(ncid, varname, varid)
    call nf95_get_var(ncid, varid, var_ana3d)
    call nf95_close(ncid)
    CALL conf_dat3d(lon_ini, lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, &
         var_ana3d)

    DO il = 1, llm_dyn
       CALL inter_barxy(lon_rad, lat_rad(:jml_dyn-1), var_ana3d(:, :, il), &
            lon_in2(:iml-1), lat_in2, var_tmp2d) 
       var_tmp3d(:, :, il) = gr_int_dyn(var_tmp2d)
    ENDDO

    ! Pour l'interpolation verticale, on interpole du haut de l'atmosphère
    ! vers le sol :
    ax = lev_dyn(llm_dyn:1:-1)
    DO ij=1, jml
       DO ii=1, iml-1
          ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
          yder = SPLINE(ax, ay)
          do il=1, lml
             var3d(ii, ij, il) = SPLINT(ax, ay, yder, pls_in(ii, ij, il))
          END do
       ENDDO
    ENDDO
    var3d(iml, :, :) = var3d(1, :, :) 

  END subroutine start_inter_3d

END MODULE startdyn
1	guez	3	MODULE startdyn
2
3			! From startvar.F, version 1.4
4			! 2006/01/27 15:14:22 Fairhead
5
6			IMPLICIT NONE
7
8			private
9			public start_init_dyn, start_inter_3d
10
11			INTEGER fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
12
13			REAL, ALLOCATABLE:: lon_ini(:), lat_ini(:)
14			! (longitude and latitude from the input file, converted to rad)
15
16			real, ALLOCATABLE:: levdyn_ini(:)
17
18			CONTAINS
19
20			SUBROUTINE start_init_dyn(tsol_2d, psol)
21
22			! Host associated variables appearing and modified in this procedure :
23			! iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn, lon_ini, lat_ini, levdyn_ini
24
25	guez	32	USE flincom, only: flininfo, flinopen_nozoom
26	guez	3	use comgeom, only: aire_2d, apoln, apols
27			use conf_dat2d_m, only: conf_dat2d
28			use inter_barxy_m, only: inter_barxy
29			use comgeom, only: rlonu, rlatv
30			use dimens_m, only: iim, jjm
31			use gr_int_dyn_m, only: gr_int_dyn
32			use start_init_orog_m, only: phis
33			use start_init_phys_m, only: start_init_phys
34	guez	39	use nr_util, only: assert, pi
35	guez	48	use netcdf, only: nf90_nowrite
36			use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid
37	guez	3
38	guez	43	REAL, intent(in):: tsol_2d(:, :) ! (iim + 1, jjm + 1)
39			REAL, intent(out):: psol(:, :) ! (iim + 1, jjm + 1) surface pressure, in Pa
40	guez	3
41			! Local:
42
43			REAL date, dt
44	guez	48	INTEGER itau(1), ncid, varid
45	guez	3	REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:)
46
47			REAL, ALLOCATABLE:: lon_dyn(:, :), lat_dyn(:, :)
48			! (longitude and latitude from the input file, in rad or degrees)
49
50	guez	36	REAL, ALLOCATABLE:: var_ana(:, :)
51			real z(iim + 1, jjm + 1)
52	guez	3	real tmp_var(iim, jjm + 1)
53
54			!--------------------------
55
56			print *, "Call sequence information: start_init_dyn"
57	guez	36	call assert((/size(tsol_2d, 1), size(psol, 1)/) == iim + 1, &
58	guez	42	"start_init_dyn size 1")
59	guez	36	call assert((/size(tsol_2d, 2), size(psol, 2)/) == jjm + 1, &
60	guez	42	"start_init_dyn size 2")
61	guez	36	CALL flininfo('ECDYN.nc', iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
62	guez	3	print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &
63			", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn
64
65			ALLOCATE(lat_dyn(iml_dyn, jml_dyn))
66			ALLOCATE(lon_dyn(iml_dyn, jml_dyn))
67			ALLOCATE(levdyn_ini(llm_dyn))
68
69	guez	32	CALL flinopen_nozoom(iml_dyn, jml_dyn, llm_dyn, &
70	guez	3	lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
71
72			ALLOCATE(var_ana(iml_dyn, jml_dyn))
73			ALLOCATE(lon_rad(iml_dyn))
74			ALLOCATE(lon_ini(iml_dyn))
75
76	guez	36	IF (MAXVAL(lon_dyn) > pi) THEN
77	guez	3	! Assume "lon_dyn" is in degrees
78	guez	36	lon_ini = lon_dyn(:, 1) * pi / 180.
79	guez	3	ELSE
80	guez	36	lon_ini = lon_dyn(:, 1)
81	guez	3	ENDIF
82
83			ALLOCATE(lat_rad(jml_dyn))
84			ALLOCATE(lat_ini(jml_dyn))
85
86	guez	36	IF (MAXVAL(lat_dyn) > pi) THEN
87			lat_ini = lat_dyn(1, :) * pi / 180.
88	guez	3	ELSE
89	guez	36	lat_ini = lat_dyn(1, :)
90	guez	3	ENDIF
91
92	guez	48	call nf95_open('ECDYN.nc', nf90_nowrite, ncid)
93
94	guez	3	! 'Z': Surface geopotential
95	guez	48	call nf95_inq_varid(ncid, 'Z', varid)
96			call nf95_get_var(ncid, varid, var_ana)
97	guez	3	CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
98			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
99			rlatv, tmp_var)
100	guez	36	z = gr_int_dyn(tmp_var)
101	guez	3
102			! 'SP': Surface pressure
103	guez	48	call nf95_inq_varid(ncid, 'SP', varid)
104			call nf95_get_var(ncid, varid, var_ana)
105	guez	3	CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
106			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
107			rlatv, tmp_var)
108	guez	36	psol = gr_int_dyn(tmp_var)
109	guez	3
110	guez	48	call nf95_close(ncid)
111
112	guez	36	psol(:iim, :) = psol(:iim, :) &
113			* (1. + (z(:iim, :) - phis(:iim, :)) / 287. / tsol_2d(:iim, :))
114	guez	3	psol(iim + 1, :) = psol(1, :)
115
116	guez	36	psol(:, 1) = SUM(aire_2d(:iim, 1) * psol(:iim, 1)) / apoln
117			psol(:, jjm + 1) = SUM(aire_2d(:iim, jjm + 1) * psol(:iim, jjm + 1)) &
118			/ apols
119	guez	3
120			END SUBROUTINE start_init_dyn
121
122			!********************************
123
124	guez	23	subroutine start_inter_3d(varname, lon_in2, lat_in2, pls_in, var3d)
125	guez	3
126	guez	36	! This procedure gets a 3D variable from a file and interpolates it.
127	guez	3
128	guez	36	use nr_util, only: assert_eq
129			use numer_rec, only: spline, splint
130	guez	3	use inter_barxy_m, only: inter_barxy
131			use gr_int_dyn_m, only: gr_int_dyn
132			use conf_dat3d_m, only: conf_dat3d
133	guez	48	use netcdf, only: nf90_nowrite
134			use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid
135	guez	3
136			CHARACTER(len=*), intent(in):: varname
137	guez	36	REAL, intent(in):: lon_in2(:) ! (iml)
138			REAL, intent(in):: lat_in2(:)
139			REAL, intent(in):: pls_in(:, :, :) ! (iml, jml, lml)
140			REAL, intent(out):: var3d(:, :, :) ! (iml, jml, lml)
141	guez	3
142			! LOCAL:
143	guez	48	INTEGER iml, jml, lml, ncid, varid
144	guez	3	INTEGER ii, ij, il
145			REAL lon_rad(iml_dyn), lat_rad(jml_dyn)
146			REAL lev_dyn(llm_dyn)
147			REAL var_tmp2d(size(lon_in2)-1, size(pls_in, 2))
148			real var_tmp3d(size(lon_in2), size(pls_in, 2), llm_dyn)
149			REAL ax(llm_dyn), ay(llm_dyn), yder(llm_dyn)
150			real var_ana3d(iml_dyn, jml_dyn, llm_dyn)
151
152			!--------------------------------
153
154			print *, "Call sequence information: start_inter_3d"
155
156	guez	23	iml = assert_eq(size(pls_in, 1), size(lon_in2), size(var3d, 1), &
157			"start_inter_3d iml")
158			jml = assert_eq(size(pls_in, 2), size(var3d, 2), "start_inter_3d jml")
159			lml = assert_eq(size(pls_in, 3), size(var3d, 3), "start_inter_3d lml")
160	guez	3
161			print *, "iml = ", iml, ", jml = ", jml
162	guez	23	print *, "varname = ", varname
163	guez	48	call nf95_open('ECDYN.nc', nf90_nowrite, ncid)
164			call nf95_inq_varid(ncid, varname, varid)
165			call nf95_get_var(ncid, varid, var_ana3d)
166			call nf95_close(ncid)
167	guez	3	CALL conf_dat3d(lon_ini, lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, &
168			var_ana3d)
169
170			DO il = 1, llm_dyn
171			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn-1), var_ana3d(:, :, il), &
172			lon_in2(:iml-1), lat_in2, var_tmp2d)
173			var_tmp3d(:, :, il) = gr_int_dyn(var_tmp2d)
174			ENDDO
175
176			! Pour l'interpolation verticale, on interpole du haut de l'atmosphère
177			! vers le sol :
178	guez	36	ax = lev_dyn(llm_dyn:1:-1)
179	guez	3	DO ij=1, jml
180			DO ii=1, iml-1
181	guez	36	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
182			yder = SPLINE(ax, ay)
183	guez	3	do il=1, lml
184	guez	24	var3d(ii, ij, il) = SPLINT(ax, ay, yder, pls_in(ii, ij, il))
185	guez	3	END do
186			ENDDO
187			ENDDO
188	guez	23	var3d(iml, :, :) = var3d(1, :, :)
189	guez	3
190	guez	23	END subroutine start_inter_3d
191	guez	3
192			END MODULE startdyn