gen_be_diags.f90

1 program gen_be_diags
2 
3    use da_control, only : stderr, stdout, filename_len
4    use da_tools_serial, only : da_get_unit
5    use da_gen_be, only : da_readwrite_be_stage2, da_readwrite_be_stage3, &
6       da_readwrite_be_stage4
7 
8    implicit none
9 
10    character*10        :: variable                   ! Variable name
11    character*8         :: uh_method                  ! Uh_method (power, scale)
12    integer             :: n_smth_sl                  ! Number of smoothing for scale-length
13    character(len=filename_len)        :: filename                   ! Input filename.
14    integer             :: nk,nk_3d                   ! Dimensions read in.
15 
16    namelist / gen_be_diags_nl / uh_method, n_smth_sl
17 
18    integer :: ounit,iunit,namelist_unit
19 
20    stderr = 0
21    stdout = 6
22 
23    call da_get_unit(ounit)
24    call da_get_unit(iunit)
25    call da_get_unit(namelist_unit)
26 
27    uh_method = 'scale'
28    n_smth_sl = 0
29 
30    open(unit=namelist_unit, file='gen_be_diags_nl.nl', &
31       form='formatted', status='old', action='read')
32    read(namelist_unit, gen_be_diags_nl)
33    close(namelist_unit)
34 
35    filename = 'be.dat'
36    open (ounit, file = filename, form='unformatted')
37 
38    write(6,'(/a)')' [1] Gather regression coefficients.'
39 
40    call da_readwrite_be_stage2( ounit, nk )
41 
42    write(6,'(/a)')' [2] Gather vertical error eigenvectors, eigenvalues.'
43 
44    variable = 'psi'
45    call da_readwrite_be_stage3( ounit, nk, variable )
46 
47    variable = 'chi_u'
48    call da_readwrite_be_stage3( ounit, nk, variable )
49 
50    variable = 't_u'
51    call da_readwrite_be_stage3( ounit, nk, variable )
52 
53    variable = 'rh'
54    call da_readwrite_be_stage3( ounit, nk, variable )
55 
56    ! To keep the dimension nk for 3d fields:
57    nk_3d = nk
58 
59    if (uh_method /= 'power') then
60       variable = 'ps_u'
61       call da_readwrite_be_stage3( ounit,  1, variable )
62    end if
63 
64    if (uh_method == 'power') then
65       write(6,'(/a)')' [3] Gather horizontal error power spectra.'
66    else if (uh_method == 'scale') then
67       write(6,'(/a)')' [3] Gather horizontal scale length.'
68    end if
69 
70    ! To assign the dimension nk for 3d fields:
71    nk = nk_3d
72 
73    variable = 'psi'
74    call da_readwrite_be_stage4( ounit, nk, uh_method, n_smth_sl, variable )
75 
76    variable = 'chi_u'
77    call da_readwrite_be_stage4( ounit, nk, uh_method, n_smth_sl, variable )
78 
79    variable = 't_u'
80    call da_readwrite_be_stage4( ounit, nk, uh_method, n_smth_sl, variable )
81 
82    variable = 'rh'
83    call da_readwrite_be_stage4( ounit, nk, uh_method, n_smth_sl, variable )
84 
85    variable = 'ps_u'
86    call da_readwrite_be_stage4( ounit, 1, uh_method, n_smth_sl, variable )
87 
88    close(ounit)
89 
90 end program gen_be_diags