/[lmdze]/trunk/dyn3d/startdyn.f
ViewVC logotype

Diff of /trunk/dyn3d/startdyn.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/startdyn.f90 revision 20 by guez, Wed Oct 15 16:19:57 2008 UTC trunk/Sources/dyn3d/startdyn.f revision 225 by guez, Mon Oct 16 12:35:41 2017 UTC
# Line 1  Line 1 
1  MODULE startdyn  MODULE startdyn
2    
3    ! From startvar.F, version 1.4    ! From startvar.F, version 1.4
4    ! 2006/01/27 15:14:22 Fairhead    ! January 27th, 2006 15:14:22
5    
6    IMPLICIT NONE    INTEGER iml_dyn, jml_dyn, llm_dyn
7    
8    private    REAL, allocatable:: lon_ini(:), lat_ini(:)
9    public start_init_dyn, start_inter_3d    ! longitude and latitude from the input file, converted to rad
10    
11    INTEGER fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn    real, allocatable:: levdyn_ini(:)
   
   REAL, ALLOCATABLE:: lon_ini(:), lat_ini(:)  
   ! (longitude and latitude from the input file, converted to rad)  
   
   real, ALLOCATABLE:: levdyn_ini(:)  
12    
13  CONTAINS  CONTAINS
14    
15    SUBROUTINE start_init_dyn(tsol_2d, psol)    SUBROUTINE start_init_dyn(tsol_2d, phis, ps)
   
     ! 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  
16    
     USE ioipsl, only: flininfo, flinopen_nozoom, flinget  
17      use comgeom, only: aire_2d, apoln, apols      use comgeom, only: aire_2d, apoln, apols
18      use conf_dat2d_m, only: conf_dat2d      use conf_dat2d_m, only: conf_dat2d
     use inter_barxy_m, only: inter_barxy  
     use comconst, only: pi  
     use comgeom, only: rlonu, rlatv  
19      use dimens_m, only: iim, jjm      use dimens_m, only: iim, jjm
20        use dynetat0_m, only: rlonu, rlatv
21      use gr_int_dyn_m, only: gr_int_dyn      use gr_int_dyn_m, only: gr_int_dyn
22      use start_init_orog_m, only: phis      use inter_barxy_m, only: inter_barxy
23      use start_init_phys_m, only: start_init_phys      use netcdf, only: nf90_nowrite
24        use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid, &
25             nf95_gw_var, find_coord
26        use nr_util, only: assert, pi
27    
28      REAL, intent(out):: tsol_2d(:, :)      REAL, intent(in):: tsol_2d(:, :) ! (iim + 1, jjm + 1)
     REAL, intent(out):: psol(:, :) ! surface pressure, in Pa  
29    
30      ! Local:      REAL, intent(in):: phis(:, :) ! (iim + 1, jjm + 1)
31        ! surface geopotential, in m2 s-2
32    
33      REAL date, dt      REAL, intent(out):: ps(:, :) ! (iim + 1, jjm + 1) surface pressure, in Pa
     INTEGER itau(1)  
     INTEGER i, j  
34    
35      CHARACTER(len=120) physfname      ! Local:
36    
37        INTEGER ncid, varid
38      REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:)      REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:)
39        REAL, ALLOCATABLE:: var_ana(:, :)
40      REAL, ALLOCATABLE:: lon_dyn(:, :), lat_dyn(:, :)      real z(iim + 1, jjm + 1)
     ! (longitude and latitude from the input file, in rad or degrees)  
   
     REAL, ALLOCATABLE:: var_ana(:, :), z(:, :)  
41      real tmp_var(iim, jjm + 1)      real tmp_var(iim, jjm + 1)
     REAL, ALLOCATABLE:: xppn(:), xpps(:)  
42    
43      !--------------------------      !--------------------------
44    
45      print *, "Call sequence information: start_init_dyn"      print *, "Call sequence information: start_init_dyn"
46      if (any((/size(tsol_2d, 1), size(psol, 1)/) /= iim + 1)) stop &      call assert((/size(tsol_2d, 1), size(phis, 1), size(ps, 1)/) == iim + 1, &
47           "start_init_phys size 1"           "start_init_dyn size 1")
48      if (any((/size(tsol_2d, 2), size(psol, 2)/) /= jjm + 1)) stop &      call assert((/size(tsol_2d, 2), size(phis, 2), size(ps, 2)/) == jjm + 1, &
49           "start_init_phys size 2"           "start_init_dyn size 2")
50      physfname = 'ECDYN.nc'  
51      print *, 'Opening the surface analysis'      call nf95_open('ECDYN.nc', nf90_nowrite, ncid)
52      CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)  
53      print *, 'Values read from "' // trim(physfname) // '":'      call find_coord(ncid, varid=varid, std_name="longitude")
54      print *, "iml_dyn = ", iml_dyn, ", jml_dyn = ", jml_dyn, &      call nf95_gw_var(ncid, varid, lon_ini)
55           ", llm_dyn = ", llm_dyn, ", ttm_dyn = ", ttm_dyn      lon_ini = lon_ini * pi / 180. ! convert to rad
56        iml_dyn = size(lon_ini)
57      ALLOCATE(lat_dyn(iml_dyn, jml_dyn))      print *, "iml_dyn = ", iml_dyn
58      ALLOCATE(lon_dyn(iml_dyn, jml_dyn))  
59      ALLOCATE(levdyn_ini(llm_dyn))      call find_coord(ncid, varid=varid, std_name="latitude")
60        call nf95_gw_var(ncid, varid, lat_ini)
61      CALL flinopen_nozoom(physfname, iml_dyn, jml_dyn, llm_dyn, &      lat_ini = lat_ini * pi / 180. ! convert to rad
62           lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)      jml_dyn = size(lat_ini)
63        print *, "jml_dyn = ", jml_dyn
64      ALLOCATE(var_ana(iml_dyn, jml_dyn))  
65      ALLOCATE(lon_rad(iml_dyn))      call nf95_inq_varid(ncid, "level", varid)
66      ALLOCATE(lon_ini(iml_dyn))      call nf95_gw_var(ncid, varid, levdyn_ini)
67        llm_dyn = size(levdyn_ini)
68      IF (MAXVAL(lon_dyn(:, :)) > pi) THEN      print *, "llm_dyn = ", llm_dyn
        ! 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  
69    
70      ALLOCATE(z(iim + 1, jjm + 1))      ALLOCATE(var_ana(iml_dyn, jml_dyn), lon_rad(iml_dyn), lat_rad(jml_dyn))
71    
72      ! 'Z': Surface geopotential      ! 'Z': Surface geopotential
73      CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)      call nf95_inq_varid(ncid, 'Z', varid)
74        call nf95_get_var(ncid, varid, var_ana)
75      CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)      CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
76      CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &      CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
77           rlatv, tmp_var)           rlatv, tmp_var)
78      z(:, :) = gr_int_dyn(tmp_var)      z = gr_int_dyn(tmp_var)
79    
80      ! 'SP': Surface pressure      ! 'SP': Surface pressure
81      CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)      call nf95_inq_varid(ncid, 'SP', varid)
82        call nf95_get_var(ncid, varid, var_ana)
83      CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)      CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, var_ana)
84      CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &      CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana, rlonu(:iim), &
85           rlatv, tmp_var)           rlatv, tmp_var)
86      psol(:, :) = gr_int_dyn(tmp_var)      ps = gr_int_dyn(tmp_var)
     CALL start_init_phys(tsol_2d)  
87    
88      ! PSOL is computed in Pascals      call nf95_close(ncid)
89    
90      DO j = 1, jjm + 1      ! Adapt the surface pressure to "phis", with the hypsometric equation:
91         DO i = 1, iim      ps(:iim, :) = ps(:iim, :) &
92            psol(i, j) = psol(i, j) &           * (1. + (z(:iim, :) - phis(:iim, :)) / 287. / tsol_2d(:iim, :))
                * (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  
93    
94      psol(:, 1) = SUM(xppn)/apoln      ps(iim + 1, :) = ps(1, :)
95      psol(:, jjm + 1) = SUM(xpps)/apols      ps(:, 1) = SUM(aire_2d(:iim, 1) * ps(:iim, 1)) / apoln
96        ps(:, jjm + 1) = SUM(aire_2d(:iim, jjm + 1) * ps(:iim, jjm + 1)) / apols
97    
98    END SUBROUTINE start_init_dyn    END SUBROUTINE start_init_dyn
99    
   !********************************  
   
   function start_inter_3d(varname, lon_in2, lat_in2, pls_in)  
   
     ! This procedure gets a 3D variable from a file and does the  
     ! interpolations needed.  
   
     USE ioipsl, 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 start_inter_3d(size(lon_in2), size(pls_in, 2), size(pls_in, 3))  
   
     ! 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), "start_inter_3d iml")  
     jml = size(pls_in, 2)  
     lml = size(pls_in, 3)  
   
     print *, "iml = ", iml, ", jml = ", jml  
     print *, "fid_dyn = ", fid_dyn, ", 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  
              start_inter_3d(ii, ij, il) &  
                   = SPLINT(ax, ay, yder, pls_in(ii, ij, il))  
           END do  
        ENDDO  
     ENDDO  
     start_inter_3d(iml, :, :) = start_inter_3d(1, :, :)  
   
   END function start_inter_3d  
   
100  END MODULE startdyn  END MODULE startdyn
Legend:
Removed from v.20  
changed lines
  Added in v.225
  ViewVC Help
Powered by ViewVC 1.1.21