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 flinget_m, only: flinget
    use comgeom, only: aire_2d, apoln, apols
    use conf_dat2d_m, only: conf_dat2d
    use inter_barxy_m, only: inter_barxy
    use comconst, only: pi
    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

    REAL, intent(out):: tsol_2d(:, :)
    REAL, intent(out):: psol(:, :) ! surface pressure, in Pa

    ! Local:

    REAL date, dt
    INTEGER itau(1)
    INTEGER i, j

    CHARACTER(len=120) physfname

    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(:, :), z(:, :)
    real tmp_var(iim, jjm + 1)
    REAL, ALLOCATABLE:: xppn(:), xpps(:)

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

    print *, "Call sequence information: start_init_dyn"
    if (any((/size(tsol_2d, 1), size(psol, 1)/) /= iim + 1)) stop &
         "start_init_phys size 1"
    if (any((/size(tsol_2d, 2), size(psol, 2)/) /= jjm + 1)) stop &
         "start_init_phys size 2"
    physfname = 'ECDYN.nc'
    print *, 'Opening the surface analysis'
    CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
    print *, 'Values read from "' // trim(physfname) // '":'
    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

    ALLOCATE(z(iim + 1, jjm + 1))

    ! 'Z': Surface geopotential
    CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, 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 flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, 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 start_init_phys(tsol_2d)

    ! PSOL is computed in Pascals

    DO j = 1, jjm + 1
       DO i = 1, iim
          psol(i, j) = psol(i, j) &
               * (1. + (z(i, j) - phis(i, j)) / 287. / tsol_2d(i, j))
       ENDDO
    ENDDO
    psol(iim + 1, :) = psol(1, :)

    ALLOCATE(xppn(iim))
    ALLOCATE(xpps(iim)) 

    DO  i = 1, iim
       xppn(i) = aire_2d( i, 1) * psol( i, 1)
       xpps(i) = aire_2d( i, jjm + 1) * psol( i, jjm + 1)
    ENDDO

    psol(:, 1) = SUM(xppn)/apoln
    psol(:, jjm + 1) = SUM(xpps)/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 does the
    ! interpolations needed.

    use flinget_m, only: flinget
    use numer_rec, only: assert_eq, 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

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

    ! LOCAL:
    INTEGER iml, jml, lml
    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
    print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &
         ", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn
    print *, 'Going into flinget to extract the 3D field.'
    CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
         var_ana3d)

    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			use flinget_m, only: flinget
27	guez	3	use comgeom, only: aire_2d, apoln, apols
28			use conf_dat2d_m, only: conf_dat2d
29			use inter_barxy_m, only: inter_barxy
30			use comconst, only: pi
31			use comgeom, only: rlonu, rlatv
32			use dimens_m, only: iim, jjm
33			use gr_int_dyn_m, only: gr_int_dyn
34			use start_init_orog_m, only: phis
35			use start_init_phys_m, only: start_init_phys
36
37			REAL, intent(out):: tsol_2d(:, :)
38			REAL, intent(out):: psol(:, :) ! surface pressure, in Pa
39
40			! Local:
41
42			REAL date, dt
43			INTEGER itau(1)
44			INTEGER i, j
45
46			CHARACTER(len=120) physfname
47
48			REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:)
49
50			REAL, ALLOCATABLE:: lon_dyn(:, :), lat_dyn(:, :)
51			! (longitude and latitude from the input file, in rad or degrees)
52
53			REAL, ALLOCATABLE:: var_ana(:, :), z(:, :)
54			real tmp_var(iim, jjm + 1)
55			REAL, ALLOCATABLE:: xppn(:), xpps(:)
56
57			!--------------------------
58
59			print *, "Call sequence information: start_init_dyn"
60			if (any((/size(tsol_2d, 1), size(psol, 1)/) /= iim + 1)) stop &
61			"start_init_phys size 1"
62			if (any((/size(tsol_2d, 2), size(psol, 2)/) /= jjm + 1)) stop &
63			"start_init_phys size 2"
64			physfname = 'ECDYN.nc'
65			print *, 'Opening the surface analysis'
66			CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
67			print *, 'Values read from "' // trim(physfname) // '":'
68			print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &
69			", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn
70
71			ALLOCATE(lat_dyn(iml_dyn, jml_dyn))
72			ALLOCATE(lon_dyn(iml_dyn, jml_dyn))
73			ALLOCATE(levdyn_ini(llm_dyn))
74
75	guez	32	CALL flinopen_nozoom(iml_dyn, jml_dyn, llm_dyn, &
76	guez	3	lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
77
78			ALLOCATE(var_ana(iml_dyn, jml_dyn))
79			ALLOCATE(lon_rad(iml_dyn))
80			ALLOCATE(lon_ini(iml_dyn))
81
82			IF (MAXVAL(lon_dyn(:, :)) > pi) THEN
83			! Assume "lon_dyn" is in degrees
84			lon_ini(:) = lon_dyn(:, 1) * pi / 180.
85			ELSE
86			lon_ini(:) = lon_dyn(:, 1)
87			ENDIF
88
89			ALLOCATE(lat_rad(jml_dyn))
90			ALLOCATE(lat_ini(jml_dyn))
91
92			IF (MAXVAL(lat_dyn(:, :)) > pi) THEN
93			lat_ini(:) = lat_dyn(1, :) * pi / 180.
94			ELSE
95			lat_ini(:) = lat_dyn(1, :)
96			ENDIF
97
98			ALLOCATE(z(iim + 1, jjm + 1))
99
100			! 'Z': Surface geopotential
101			CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
102			CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
103			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
104			rlatv, tmp_var)
105			z(:, :) = gr_int_dyn(tmp_var)
106
107			! 'SP': Surface pressure
108			CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
109			CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
110			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
111			rlatv, tmp_var)
112			psol(:, :) = gr_int_dyn(tmp_var)
113			CALL start_init_phys(tsol_2d)
114
115			! PSOL is computed in Pascals
116
117			DO j = 1, jjm + 1
118			DO i = 1, iim
119			psol(i, j) = psol(i, j) &
120			* (1. + (z(i, j) - phis(i, j)) / 287. / tsol_2d(i, j))
121			ENDDO
122			ENDDO
123			psol(iim + 1, :) = psol(1, :)
124
125			ALLOCATE(xppn(iim))
126			ALLOCATE(xpps(iim))
127
128			DO i = 1, iim
129			xppn(i) = aire_2d( i, 1) * psol( i, 1)
130			xpps(i) = aire_2d( i, jjm + 1) * psol( i, jjm + 1)
131			ENDDO
132
133			psol(:, 1) = SUM(xppn)/apoln
134			psol(:, jjm + 1) = SUM(xpps)/apols
135
136			END SUBROUTINE start_init_dyn
137
138			!********************************
139
140	guez	23	subroutine start_inter_3d(varname, lon_in2, lat_in2, pls_in, var3d)
141	guez	3
142			! This procedure gets a 3D variable from a file and does the
143			! interpolations needed.
144
145	guez	32	use flinget_m, only: flinget
146	guez	13	use numer_rec, only: assert_eq, spline, splint
147	guez	3	use inter_barxy_m, only: inter_barxy
148			use gr_int_dyn_m, only: gr_int_dyn
149			use conf_dat3d_m, only: conf_dat3d
150
151			CHARACTER(len=*), intent(in):: varname
152			REAL, intent(in):: lon_in2(:), lat_in2(:)
153			REAL, intent(in):: pls_in(:, :, :)
154	guez	23	REAL, intent(out):: var3d(:, :, :)
155	guez	3
156			! LOCAL:
157			INTEGER iml, jml, lml
158			INTEGER ii, ij, il
159			REAL lon_rad(iml_dyn), lat_rad(jml_dyn)
160			REAL lev_dyn(llm_dyn)
161			REAL var_tmp2d(size(lon_in2)-1, size(pls_in, 2))
162			real var_tmp3d(size(lon_in2), size(pls_in, 2), llm_dyn)
163			REAL ax(llm_dyn), ay(llm_dyn), yder(llm_dyn)
164			real var_ana3d(iml_dyn, jml_dyn, llm_dyn)
165
166			!--------------------------------
167
168			print *, "Call sequence information: start_inter_3d"
169
170	guez	23	iml = assert_eq(size(pls_in, 1), size(lon_in2), size(var3d, 1), &
171			"start_inter_3d iml")
172			jml = assert_eq(size(pls_in, 2), size(var3d, 2), "start_inter_3d jml")
173			lml = assert_eq(size(pls_in, 3), size(var3d, 3), "start_inter_3d lml")
174	guez	3
175			print *, "iml = ", iml, ", jml = ", jml
176	guez	23	print *, "varname = ", varname
177	guez	3	print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &
178			", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn
179	guez	20	print *, 'Going into flinget to extract the 3D field.'
180	guez	3	CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
181			var_ana3d)
182
183			CALL conf_dat3d(lon_ini, lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, &
184			var_ana3d)
185
186			DO il = 1, llm_dyn
187			CALL inter_barxy(lon_rad, lat_rad(:jml_dyn-1), var_ana3d(:, :, il), &
188			lon_in2(:iml-1), lat_in2, var_tmp2d)
189			var_tmp3d(:, :, il) = gr_int_dyn(var_tmp2d)
190			ENDDO
191
192			! Pour l'interpolation verticale, on interpole du haut de l'atmosphère
193			! vers le sol :
194			ax(:) = lev_dyn(llm_dyn:1:-1)
195			DO ij=1, jml
196			DO ii=1, iml-1
197			ay(:) = var_tmp3d(ii, ij, llm_dyn:1:-1)
198			yder(:) = SPLINE(ax, ay)
199			do il=1, lml
200	guez	24	var3d(ii, ij, il) = SPLINT(ax, ay, yder, pls_in(ii, ij, il))
201	guez	3	END do
202			ENDDO
203			ENDDO
204	guez	23	var3d(iml, :, :) = var3d(1, :, :)
205	guez	3
206	guez	23	END subroutine start_inter_3d
207	guez	3
208			END MODULE startdyn