text:
module mod_test_parameters
use blas_sparse_namedconstants
use properties
! *** The test matrices
      integer, parameter :: MATRIX_A = 1
      integer, parameter :: MATRIX_A_SU = 2
      integer, parameter :: MATRIX_A_SL = 3
      integer, parameter :: MATRIX_T = 4
      integer, parameter :: MATRIX_U = 5
      integer, parameter :: MATRIX_B = 6
      integer, parameter :: MATRIX_X = 7
      integer, parameter :: MATRIX_Y = 8
      integer, parameter :: MATRIX_Y_SU = 9
      integer, parameter :: MATRIX_Y_SL = 10
      integer, parameter :: MATRIX_M = 11
      integer, parameter :: MATRIX_N = 12
      integer, parameter :: MATRIX_Z = 13
      integer, parameter :: MATRIX_C = 14
      integer, parameter :: MATRIX_I = 15
      integer, parameter :: MATRIX_J = 16
! *** Which functionality of MV/MM mult./Tri solve should be called ?
      integer, parameter :: O_MATRIX = 0 !original
      integer, parameter :: T_MATRIX = 1 !transpose
      integer, parameter :: H_MATRIX = 2 !hermit

! *** Level 1 routines
real(kind=dp),dimension(3)::x_dot=(/1.,3.,5./)
real(kind=dp),dimension(6)::y_dot=(/2.,4.,6.,8.,10.,12./)
integer ,dimension(3),parameter ::indx_dot=(/2,3,4/)

real(kind=dp),dimension(3)::x_axpy=(/51.,53.,55./)
real(kind=dp),dimension(6)::y_axpy=(/52.,54.,56.,58.,60.,62./)
integer ,dimension(3)::indx_axpy=(/2,3,4/)
real(kind=dp),parameter::alpha=2.

real(kind=dp),dimension(3)::x_ga=(/31.,33.,35./)
real(kind=dp),dimension(6)::y_ga=(/32.,34.,36.,38.,40.,42./)
integer ,dimension(3)::indx_ga=(/2,3,4/)

real(kind=dp),dimension(3)::x_gz=(/11.,13.,15./)
real(kind=dp),dimension(6)::y_gz=(/12.,14.,16.,18.,20.,22./)
integer ,dimension(3) ::indx_gz=(/2,3,4/)

real(kind=dp),dimension(3)::x_sc=(/41.,43.,45./)
real(kind=dp),dimension(6)::y_sc=(/42.,44.,46.,48.,50.,52./)
integer ,dimension(3) ::indx_sc=(/2,3,4/)

complex(kind=dp),dimension(3)::xc_dot=(/(1.,1.),(3.,1.),(5.,1.)/)
complex(kind=dp),dimension(6)::yc_dot=(/(2.,1.),(4.,1.),(6.,1.),(8.,1.),(10.,1.),(12.,1.)/)
integer ,dimension(3),parameter ::indxc_dot=(/2,3,4/)

complex(kind=dp),dimension(3)::xc_axpy=(/(51.,1.),(53.,1.),(55.,1.)/)
complex(kind=dp),dimension(6)::yc_axpy=(/(52.,1.),(54.,1.),(56.,1.),(58.,1.),(60.,1.),(62.,1.)/)
integer ,dimension(3)::indxc_axpy=(/2,3,4/)
complex(kind=dp),parameter::alphac=(2.,1)

complex(kind=dp),dimension(3)::xc_ga=(/(31.,1.),(33.,1.),(35.,1.)/)
complex(kind=dp),dimension(6)::yc_ga=(/(32.,1.),(34.,1.),(36.,1.),(38.,1.),(40.,1.),(42.,1.)/)
integer ,dimension(3)::indxc_ga=(/2,3,4/)

complex(kind=dp),dimension(3)::xc_gz=(/(11.,1.),(13.,1.),(15.,1.)/)
complex(kind=dp),dimension(6)::yc_gz=(/(12.,1.),(14.,1.),(16.,1.),(18.,1.),(20.,1.),(22.,1.)/)
integer ,dimension(3) ::indxc_gz=(/2,3,4/)

complex(kind=dp),dimension(3)::xc_sc=(/(41.,1.),(43.,1.),(45.,1.)/)
complex(kind=dp),dimension(6)::yc_sc=(/(42.,1.),(44.,1.),(46.,1.),(48.,1.),(50.,1.),(52.,1.)/)
integer ,dimension(3) ::indxc_sc=(/2,3,4/)

! *** Level 2,3 and handle management routines
 !-----------------------------------------------------------------
      ! Test - Matrices :
      !
      !    /11  0 13 14  0\     /11  0 13 14  0\     /11  0 31 41 51\
      !    | 0  0 23 24  0|     | 0  0 23 24  0|     | 0  0 32 42 52|
      ! A= |31 32 33 34  0|A_SU=|13 23 33 34  0|A_SL=|31 32 33  0  0|
      !    | 0 42  0 44  0|     |14 24  0 44  0|     |41 42  0 44  0|
      !    \51 52  0  0 55/     \ 0  0  0  0 55/     \51 52  0  0 55/
      !
      !    / 1  4  7 10  13\
      ! B= | 2  5  8 11  14|
      !    \ 3  6  9 12  15/
      !
      !    / 1 i\
      ! C= |    |
      !    \-i 1/
      !
      !    / 1 i\
      ! D= |    |
      !    \-i 1/
      !
      ! I = A (integer storage)
      !
      ! J = T (integer storage)
      !
      !    /  1  2 |  4  7 \     /  1  0 |      \
      !    |  0  3 |  5  8 |     |  2  5 |      |
      ! M= |-------+-------|  N= |-------+------|
      !    |       |  6  9 |     |  3  6 | 8  0 |
      !    \       |  0 10 /     \  4  7 | 9 10 /
      !
      !    / 1  1  1  1  1\     / 1            \
      !    |    1  1  1  1|     | 1  1         |
      ! T= |       1  1  1|   U=| 1  1  1      |
      !    |          1   |     | 1  1  1  1   |
      !    \             1/     \ 1  1  1     1/
      !
      !    /11 12 |  0  0 | 15 16|  0  0\
      !    |21 22 |  0  0 | 25 26|  0  0|
      !    |------+-------+------+------|
      ! X= | 0  0 | 33  0 | 35 36|  0  0|
      !    | 0  0 | 43 44 | 45 46|  0  0|
      !    |------+-------+------+------|
      !    |51 52 |  0  0 |  0  0|  0  0|
      !    \61 62 |  0  0 |  0  0|  0  0/
      !
      !    / 1  3 |  5  7 |  0  0\
      !    | 3  4 |  6  8 |  0  0|
      !    |------+-------+------|
      ! Y= | 9 11 | 13 15 |  0  0|
      !    |10 12 | 15 16 |  0  0|
      !    |------+-------+------|
      !    | 0  0 |  0  0 |  0  0|
      !    \ 0  0 |  0  0 |  0  0/
      !
      !    / 4  2 |  0  0  0 |  1 |  0  0  0 | -1  1 \
      !    | 1  5 |  0  0  0 |  2 |  0  0  0 |  0 -1 |
      !    |------+----------+----+----------+-------|
      !    | 0  0 |  6  1  2 |  2 |  0  0  0 |  0  0 |
      !    | 1  5 |  2  7  1 |  0 |  0  0  0 |  0  0 |
      !    | 1  5 | -1  2  9 |  3 |  0  0  0 |  0  0 |
      !    |------+----------+----+----------+-------|
      ! Z= | 2  1 |  3  4  5 | 10 |  4  3  2 |  0  0 |
      !    |------+----------+----+----------+-------|
      !    | 0  0 |  0  0  0 |  4 | 13  4  2 |  0  0 |
      !    | 0  0 |  0  0  0 |  3 |  3 11  3 |  0  0 |
      !    | 0  0 |  0  0  0 |  0 |  2  0  7 |  0  0 |
      !    |------+----------+----+----------+-------|
      !    | 8  4 |  0  0  0 |  0 |  0  0  0 | 25  3 |
      !    \-2  3 |  0  0  0 |  0 |  0  0  0 |  8 12 /
      !
      !-----------------------------------------------------------------

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      !    /11  0 13 14  0\
      !    | 0  0 23 24  0|
      ! A= |31 32 33 34  0|
      !    | 0 42  0 44  0|
      !    \51 52  0  0 55/
      !
   real(kind=dp), dimension(14)  ::&
 A_val=(/ 11.,51.,31.,32.,34.,52.,13.,23.,33.,14.,24.,42.,55.,44. /)
      integer, dimension(14)  ::&
 A_indx= (/ 1,5,3,3,3,5,1,2,3,1,2,4,5,4 /)
      integer, dimension(14)  ::&
 A_jndx = (/ 1,1,1,2,4,2,3,3,3,4,4,2,5,4 /)
     integer ::A_m=5
     integer ::A_n=5

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      !      /11  0 31  0 51\
      !      | 0  0 32 42 52|
      ! A_SL=|31 32 33  0  0|
      !      | 0 42  0 44  0|
      !      \51 52  0  0 55/
      !

   real(kind=dp),dimension(5) ::A_SL_1=(/ 11., 0.,31., 0.,51./)
   integer,dimension(5) ::A_SL_1_INDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::A_SL_2=(/0.,  0.,32., 42., 52./)
   integer,dimension(5) ::A_SL_2_INDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::A_SL_3=(/31.,32.,33., 0.,0./)
   integer,dimension(5) ::A_SL_3_INDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::A_SL_4=(/0.,42.,0., 44.,0./)
   integer,dimension(5) ::A_SL_4_INDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::A_SL_5=(/51.,52.,0., 0.,55./)
   integer,dimension(5) ::A_SL_5_INDX=(/1,2,3,4,5/)

   integer ::A_SL_m=5
   integer,parameter::A_SL_n=5

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   !     /11  0 13 14  0\
   !     | 0  0 23 24  0|
   !A_SU=|13 23 33 34  0|
   !     |14 24  0 44  0|
   !     \ 0  0  0  0 55/

      real(kind=dp),dimension(5,5),parameter::&
           A_SU=reshape((/11.,0.,0.,0.,0.,0.,0.,0.,0.,0.,13.,23.,33.,0.,0.,14.,24.,34.,44.,0.,0.,0.,0.,0.,55./),shape=(/5,5/))
      integer,dimension(5) ::A_SU_indx=(/1,2,3,4,5/)
      integer,dimension(5) ::A_SU_jndx=(/1,2,3,4,5/)
      integer,parameter::A_SU_m=5
      integer,parameter::A_SU_n=5

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      !    / 1  4  7 10  13\
      ! B= | 2  5  8 11  14|
      !    \ 3  6  9 12  15/
      !

   real(kind=dp),dimension(5) ::B_1=(/ 1.,4.,7.,10.,13./)
   integer,dimension(5) ::B_1_JNDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::B_2=(/2.,5.,8.,11.,14./)
   integer,dimension(5) ::B_2_JNDX=(/1,2,3,4,5/)

   real(kind=dp),dimension(5) ::B_3=(/3.,  6.,9.,12.,15./)
   integer,dimension(5) ::B_3_JNDX=(/1,2,3,4,5/)
   integer ::B_m=3
   integer ::B_n=5

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1 0 i\
      !C= |  0 0 0 |
      !    \-i 0 1/

complex(kind=dp),dimension(4),parameter ::VAL_C=(/(1.,0.),(0.,1.),(0.,-1.),(1.,0.)/)
integer ,dimension(4)  ::INDX_C=(/1,1,3,3/)
integer ,dimension(4)  ::JNDX_C=(/1,3,1,3/)
integer,parameter::C_m=3
integer,parameter::C_n=3

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1 -i -i\
      !D= |  i  1 -i |
      !    \ i  i  1 /

complex(kind=dp),dimension(9),parameter ::VAL_d=(/(1.,0.),(0.,-1.),(0.,-1.),(0.,1.),(1.,0.),(0.,-1.),(0.,1.),(0.,1.),(1.,0.)/)
integer ,dimension(9)  ::INDX_d=(/1,1,1,2,2,2,3,3,3/)
integer ,dimension(9)  ::JNDX_d=(/1,2,3,1,2,3,1,2,3/)
integer,parameter::d_m=3
integer,parameter::d_n=3

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1 -i -i \
     !DT= |     1  -i |
      !    \       1 /

complex(kind=dp),dimension(6),parameter ::VAL_DT=(/(1.,0.),(0.,-1.),(0.,-1.),(1.,0.),(0.,-1.),(1.,0.)/)
integer ,dimension(6)  ::INDX_DT=(/1,1,1,2,2,3/)
integer ,dimension(6)  ::JNDX_DT=(/1,2,3,2,3,3/)
integer,parameter::DT_m=3
integer,parameter::DT_N=3
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    /11  0 13 14  0\
      !    | 0  0 23 24  0|
      ! I= |31 32 33 34  0|
      !    | 0 42  0 44  0|
      !    \51 52  0  0 55/
      !
   integer, dimension(14)  ::&
 I_val=(/ 11,51,31,32,34,52,13,23,33,14,24,42,55,44 /)
      integer, dimension(14)  ::&
 I_indx= (/ 1,5,3,3,3,5,1,2,3,1,2,4,5,4 /)
      integer, dimension(14)  ::&
 I_jndx = (/ 1,1,1,2,4,2,3,3,3,4,4,2,5,4 /)
     integer ::I_m=5
     integer ::I_n=5

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1  1  1  1  1\
      !    |    1  1  1  1|
      ! J= |       1  1  1|
      !    |          1   |
      !    \             1/
      !

integer,dimension(14) ::J_VAL=1
integer,dimension(14) ::J_indx=(/1,1,2,1,2,3,1,2,3,4,1,2,3,5/)
integer,dimension(14) ::J_jndx=(/1,2,2,3,3,3,4,4,4,4,5,5,5,5/)
integer,parameter::J_m=5
integer,parameter::J_n=5
integer,parameter::J_nz=14

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    /  1  2 |  4  7 \
      !    |  0  3 |  5  8 |
      ! M= |-------+-------|
      !    |       |  6  9 |
      !    \       |  0 10 /

real(kind=dp),dimension(2,2), parameter::&
 M11=reshape((/1,0,2,3/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 M12=reshape((/4,5,7,8/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
M22=reshape((/6,0,9,10/),shape=(/2,2/))

integer,parameter::M_Mb=2
integer,parameter::M_Nb=2
integer,parameter::M_k=2
integer,parameter::M_l=2
      !
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

!    /  1  0 |      \
!    |  2  5 |      |
! N= |-------+------|
!    |  3  6 | 8  0 |
!    \  4  7 | 9 10 /

real(kind=dp),dimension(2,2), parameter::&
 N11=reshape((/1,2,0,5/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 N21=reshape((/3,4,6,7/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
N22=reshape((/8,9,0,10/),shape=(/2,2/))

integer,parameter::N_Mb=2
integer,parameter::N_Nb=2
integer,parameter::N_k=2
integer,parameter::N_l=2
      !
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1  1  1  1  1\
      !    |    1  1  1  1|
      ! T= |       1  1  1|
      !    |          1   |
      !    \             1/
      !

real(kind=dp),dimension(14) ::T_VAL=1.

integer,dimension(14) ::T_indx=(/1,1,2,1,2,3,1,2,3,4,1,2,3,5/)
integer,dimension(14) ::T_jndx=(/1,2,2,3,3,3,4,4,4,4,5,5,5,5/)
integer,parameter::T_m=5
integer,parameter::T_n=5
integer,parameter::T_nz=14

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

!    / 1            \
!    | 1  1         |
!  U=| 1  1  1      |
!    | 1  1  1  1   |
!    \ 1  1  1     1/

real(kind=dp),dimension(14) ::U_VAL=1.

integer,dimension(14) ::U_jndx=(/1,1,2,1,2,3,1,2,3,4,1,2,3,5/)
integer,dimension(14) ::U_indx=(/1,2,2,3,3,3,4,4,4,4,5,5,5,5/)
integer,parameter::U_m=5
integer,parameter::U_n=5
integer,parameter::U_nz=14

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    /11 12 |  0  0 | 15 16|  0  0\
      !    |21 22 |  0  0 | 25 26|  0  0|
      !    |------+-------+------+------|
      ! X= | 0  0 | 33  0 | 35 36|  0  0|
      !    | 0  0 | 43 44 | 45 46|  0  0|
      !    |------+-------+------+------|
      !    |51 52 |  0  0 |  0  0|  0  0|
      !    \61 62 |  0  0 |  0  0|  0  0/

real(kind=dp),dimension(2,2), parameter::&
 X11=reshape((/11,21,12,22/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 X31=reshape((/51,61,52,62/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
X22=reshape((/33,43,0,44/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
X13=reshape((/15,25,16,26/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
X23=reshape((/35,45,36,46/),shape=(/2,2/))
integer,parameter::X_Mb=3
integer,parameter::X_Nb=4
integer,parameter::X_k=2
integer,parameter::X_l=2

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1  3 |  5  7 |  0  0\
      !    | 3  4 |  6  8 |  0  0|
      !    |------+-------+------|
      ! Y= | 9 11 | 13 15 |  0  0|
      !    |10 12 | 15 16 |  0  0|
      !    |------+-------+------|
      !    | 0  0 |  0  0 |  0  0|
      !    \ 0  0 |  0  0 |  0  0/

real(kind=dp),dimension(2,2), parameter::&
 Y11=reshape((/1,3,3,4/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 Y21=reshape((/9,10,11,12/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y12=reshape((/5,6,7,8/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y22=reshape((/13,15,15,16/),shape=(/2,2/))

integer,parameter::Y_Mb=3
integer,parameter::Y_Nb=3
integer,parameter::Y_k=2
integer,parameter::Y_l=2

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !       / 1  3 |  5  7 |  0  0\
      !       | 3  4 |  6  8 |  0  0|
      !       |------+-------+------|
      ! Y_SU= | 5  6 | 13 15 |  0  0|
      !       | 7  8 | 15 16 |  0  0|
      !       |------+-------+------|
      !       | 0  0 |  0  0 |  0  0|
      !       \ 0  0 |  0  0 |  0  0/

real(kind=dp),dimension(2,2), parameter::&
 Y_SU11=reshape((/1,3,3,4/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 Y_SU21=reshape((/5,7,6,8/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y_SU12=reshape((/5,6,7,8/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y_SU22=reshape((/13,15,15,16/),shape=(/2,2/))

integer,parameter::Y_SU_Mb=3
integer,parameter::Y_SU_Nb=3
integer,parameter::Y_SU_k=2
integer,parameter::Y_SU_l=2

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !    / 1  3 |  9 10 |  0  0\
      !    | 3  4 | 11 12 |  0  0|
      !    |------+-------+------|
    !Y_SL= | 9 11 | 13 15 |  0  0|
      !    |10 12 | 15 16 |  0  0|
      !    |------+-------+------|
      !    | 0  0 |  0  0 |  0  0|
      !    \ 0  0 |  0  0 |  0  0/

real(kind=dp),dimension(2,2), parameter::&
 Y_SL11=reshape((/1,3,3,4/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter::&
 Y_SL21=reshape((/9,10,11,12/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y_SL12=reshape((/9,11,10,12/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
Y_SL22=reshape((/13,15,15,16/),shape=(/2,2/))

integer,parameter::Y_SL_Mb=3
integer,parameter::Y_SL_Nb=3
integer,parameter::Y_SL_k=2
integer,parameter::Y_SL_l=2

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      !  test for Variable blocks
      !
      !    / 4  2 |  0  0  0 |  1 |  0  0  0 | -1  1 \
      !    | 1  5 |  0  0  0 |  2 |  0  0  0 |  0 -1 |
      !    |------+----------+----+----------+-------|
      !    | 0  0 |  6  1  2 |  2 |  0  0  0 |  0  0 |
      !    | 1  5 |  2  7  1 |  0 |  0  0  0 |  0  0 |
      !    | 1  5 | -1  2  9 |  3 |  0  0  0 |  0  0 |
      !    |------+----------+----+----------+-------|
      ! Z= | 2  1 |  3  4  5 | 10 |  4  3  2 |  0  0 |
      !    |------+----------+----+----------+-------|
      !    | 0  0 |  0  0  0 |  4 | 13  4  2 |  0  0 |
      !    | 0  0 |  0  0  0 |  3 |  3 11  3 |  0  0 |
      !    | 0  0 |  0  0  0 |  0 |  2  0  7 |  0  0 |
      !    |------+----------+----+----------+-------|
      !    | 8  4 |  0  0  0 |  0 |  0  0  0 | 25  3 |
      !    \-2  3 |  0  0  0 |  0 |  0  0  0 |  8 12 /
      !
      !-----------------------------------------------------------------
real(kind=dp),dimension(2,2), parameter::&
 Z11=reshape((/4.0,1.0,2.0,5.0/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
 Z15=reshape((/-1.0,0.0,1.0,-1.0/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
 Z55=reshape((/25.0,8.0,3.0,12.0/),shape=(/2,2/))
real(kind=dp),dimension(2,2), parameter:: &
 Z51=reshape((/8.0,-2.0,4.0,3.0/),shape=(/2,2/))
real(kind=dp),dimension(3,3), parameter:: &
Z22=reshape((/6.0,2.0,-1.0,1.0,7.0,2.0,2.0,1.0,9.0/),shape=(/3,3/))
real(kind=dp),dimension(3,3), parameter:: &
Z44=reshape((/13.0,3.0,2.0,4.0,11.0,0.0,2.0,3.0,7.0/),shape=(/3,3/))
real(kind=dp),dimension(2,1), parameter:: &
Z13=reshape((/1.0,2.0/),shape=(/2,1/))
real(kind=dp),dimension(1,3), parameter ::&
Z32=reshape((/3.0,4.0,5.0/),shape=(/1,3/))
real(kind=dp),dimension(1,3), parameter ::&
Z34=reshape((/4.0,3.0,2.0/),shape=(/1,3/))
real(kind=dp),dimension(3,1), parameter ::&
Z43=reshape((/4.0,3.0,0.0/),shape=(/3,1/))
real(kind=dp),dimension(3,1), parameter ::&
Z23=reshape((/2.0,0.0,3.0/),shape=(/3,1/))
real(kind=dp),dimension(1,2), parameter ::&
Z31=reshape((/2.0,1.0/),shape=(/1,2/))
real(kind=dp),dimension(1,1), parameter ::&
 Z33=reshape((/10.0/),shape=(/1,1/))

integer ,dimension(5)  ::Z_kk=(/2,3,1,3,2/)
integer ,dimension(5)  ::Z_ll=(/2,3,1,3,2/)
integer,parameter::Z_Mb=5
integer,parameter::Z_Nb=5

     !---------------------------------------------------!
      ! Test - Results for multiplication with x(i)=dbl(i)!
      !---------------------------------------------------!

      real(kind=dp) :: res_usdot = 62.
      real(kind=dp) :: res_usaxpy(6) =(/52.,156.,162.,168.,60.,62./)
      real(kind=dp) :: res_usga(3) =(/34.,36.,38./)
      real(kind=dp) :: res_usgz_x(3)=(/14.,16.,18./)
      real(kind=dp) :: res_usgz_y(6)=(/12.,0.,0.,0.,20.,22./)
      real(kind=dp) :: res_ussc(6)=(/42.,41.,43.,45.,50.,52./)
      complex(kind=dp)  :: resc_usdot=(65.,-9.)
      complex(kind=dp) :: resc_usaxpy(6) =(/(52.,1.),(155.,54.),(161.,56.),(167.,58.),(60.,1.),(62.,1.)/)
      complex(kind=dp) :: resc_usga(3) =(/(34.,1.),(36.,1.),(38.,1.)/)
      complex(kind=dp) :: resc_usgz_x(3)=(/(14.,1.),(16.,1.),(18.,1.)/)
      complex(kind=dp) :: resc_usgz_y(6)=(/(12.,1.),(0.,0.),(0.,0.),(0.,0.),(20.,1.),(22.,1.)/)
      complex(kind=dp) :: resc_ussc(6)=(/(42.,1.),(41.,1.),(43.,1.),(45.,1.),(50.,1.),(52.,1.)/)
      real(kind=dp) :: Ax(5) =(/106.,165.,330.,260.,430./)
      real(kind=dp) :: A_SUx(5)= (/106.,165.,294.,340.,275./)
      real(kind=dp) :: A_SLx(5) = (/359.,524.,194.,260.,430./)
      real(kind=dp) :: ATx(5) = (/359.,524.,158.,340.,275./)
      real(kind=dp) :: AT_SLx(5) = (/106.,165.,294.,340.,275./)
      real(kind=dp) :: AT_SUx(5) = (/359.,524.,194.,260.,430./)
      real(kind=dp) :: Bx(3) = (/ 135.,150.,165. /)
      real(kind=dp) :: BTx(5) = (/14.,32.,50.,68.,86./)
      real(kind=dp) :: Mx(4) = (/ 45.,53.,54.,40./)
      real(kind=dp) :: MTx(4) = (/1.,8.,32.,90./)
      real(kind=dp) :: Nx(4) = (/ 1.,12.,39.,85./)
      real(kind=dp) :: NTx(4) = (/30.,56.,60.,40./)
      real(kind=dp) :: Tx(5) =(/15.,14.,12.,4.,5./)
      real(kind=dp) :: TTx(5) = (/ 1.,3.,6.,10.,11. /)
      real(kind=dp) :: Ux(5) = (/ 1.,3.,6.,10.,11. /)
      real(kind=dp) :: UTx(5) = (/15.,14.,12.,4.,5./)
      real(kind=dp) :: Xx(6) =  (/206.,346.,490.,806.,155.,185./)
      real(kind=dp) :: XTx(8) = (/674.,688.,271.,176.,350.,360.,0.,0./)
      real(kind=dp) :: Yx(6) =  (/50.,61.,130.,143.,0.,0./)
      real(kind=dp) :: YTx(6) = (/74.,92.,116.,132.,0.,0./)
      real(kind=dp) :: Y_SUx(6) = (/50.,61.,116.,132.,0.,0./)
      real(kind=dp) :: Y_SLx(6) = (/74.,92.,130.,143.,0.,0./)
      real(kind=dp) :: Zx(11) = (/15.,12.,44.,39.,68.,184.,165.,154.,77.,299.,216./)
      real(kind=dp) :: ZTx(11) = (/76.,91.,39.,65.,85.,138.,157.,134.,113.,337.,161./)
      complex(kind=dp) :: Cx(3) =(/(0.,2.),(0.,0.),(2.,0.)/)
      complex(kind=dp) :: CTx(3)=(/(0.,2.),(0.,0.),(2.,0.)/)
      complex(kind=dp) :: CHx(3)=(/(0.,2.),(0.,0.),(2.,0.)/)
      complex(kind=dp) :: dx(3) =(/(3.,-1.),(1.,1.),(-1.,3.)/)
      complex(kind=dp) :: dTx(3)=(/(0.,2.),(0.,0.),(2.,0.)/)
      complex(kind=dp) :: dHx(3)=(/(0.,2.),(0.,0.),(2.,0.)/)

      complex(kind=dp) :: DTTx(3) =(/(3.,-1.),(2.,0.),(1.,1.)/)
      complex(kind=dp) :: DTHx(3) =(/(3.,-1.),(1.,1.),(-1.,3.)/)

      integer :: Ix(5) =(/106,165,330,260,430/)
      integer :: Jx(5) =(/15,14,12,4,5/)

end module mod_test_parameters

max_stars_repo_name: ittnas/qsc

max_issues_repo_name: ittnas/qsc

max_forks_repo_name: ittnas/qsc


text:
c--- This subroutine is a basic implementaion of hep-ph/9801442 valid for
c---- f(p1)+f(p2) --> gamma(p3)+f(p4)+f(p5) it should only be used for checking the more general
c---- routine frixione.f C Williams 9th Nov 2010

      subroutine basic_frix(pin,passed,isub)
      implicit none
      include 'constants.f'
      include 'frag.f'
      include 'process.f'
      include 'part.f'
      double precision pin(mxpart,4)
      integer isub
      logical passed,in_cone
      double precision ret_ET,R
      double precision R34,R35,ET3,ET4,ET5,alpha

      passed=.true.

c--- This routine also called for 'gamjet' -> always need second branch
      if ( ((part .eq. 'lord') .or. (part .eq. 'virt')
     &.or. (isub .eq. 1)) .and. (case .eq. 'dirgam') ) then
c---  LO, virtual and dipole pieces for 'dirgam'
         ET3=ret_ET(pin,3)
         ET4=ret_ET(pin,4)

         R34=R(pin,3,4)

         alpha=ET3/(1d0-dcos(cone_ang))

c---- Calculate prefac for check
         if (R34 .lt. cone_ang) then
           passed=in_cone(R34,ET4,alpha)
         else
           passed = .true.
         endif
         return

c---   isub=0 real piece for 'dirgam', always for 'gamjet'
      elseif (isub .eq. 0) then

c---- Calculate angles and energies
         ET3=ret_ET(pin,3)
         ET4=ret_ET(pin,4)
         ET5=ret_ET(pin,5)

         R34=R(pin,3,4)
         R35=R(pin,3,5)

         alpha=ET3/(1d0-dcos(cone_ang))

         if((R34 .lt. cone_ang) .and. (R35 .lt. cone_ang)) then
c----     2 partons cannot be unresolved at NLO
            passed = .false.
c---- This code is not needed
c           if (R34 .lt. R35) then
c             minR=R34
c             maxR=R35
c             minET=ET4
c           else
c             minR=R35
c             maxR=R34
c             minET=ET5
c           endif
c           passed=in_cone(minR,minET,alpha)
c           if (passed) passed=in_cone(maxR,ET4+ET5,alpha)
         elseif (R34 .lt. cone_ang) then
            passed=in_cone(R34,ET4,alpha)
         elseif (R35 .lt. cone_ang) then
            passed=in_cone(R35,ET5,alpha)
         else
            passed = .true.
         endif

      endif

      return
      end

      logical function in_cone(Rij,Ejet,pref)
      implicit none
      double precision Rij,Ejet,pref

      if ( Ejet .lt. pref*(1d0-dcos(Rij))) then
         in_cone = .true.
      else
         in_cone = .false.
      endif

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
interface good7
   subroutine s7a(x,y,z)
     class(pear) :: x,y,z
   end subroutine s7a
   subroutine s7b(x,z,w)
     class(fruit) :: x
     class(bosc) :: z
     class(apple), optional :: w
   end subroutine s7b
end interface good7

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
c               Copyright 1993 Colorado State University
c                       All Rights Reserved

      subroutine csched(cflow,protop,probot,
     &                  froma1,tob1,
     &                  froma2,tob2,
     &                  frac,accum)

      implicit none
      include 'const.inc'
      include 'zztim.inc'

c ... Argument declarations
      double precision cflow, protop, probot, froma1, tob1, froma2,
     &        tob2, frac,accum(2)

c ... Schedule C flows for decomposition from Box A to Box B
c ... written by vek 05/91

c ... Input:
c ...   cflow  = total C flow from Box A to Box B
c ...   protop = proportion top, either the labeled pool
c ...            OR cisofr or cisotf
c ...   probot = proportion bottom, either the total pool or 1
c ...   frac   = amount of fractionation
c ...   time   = simulation time (passed in /zztim/)
c ...
c ... Transput:
c ...   froma1   = state variable for unlabeled C in Box A
c ...   tob1     = state variable for unlabeled C in Box B
c ...   froma2   = state variable for labeled C in Box A
c ...   tob2     = state variable for labeled C in Box B
c ...   accum(1) = unlabeled accumulator
c ...   accum(2) = labeled accumulator

c ... Local variables
      double precision cflow1, cflow2

c ... Determine the amount of labeled to flow
      cflow2 = cflow * (protop / probot) * frac

c ... Determine the amount of unlabeled to flow
      cflow1 = cflow - cflow2

c ... Flow the amounts
      call flow(froma1,tob1,time,cflow1)
      call flow(froma2,tob2,time,cflow2)

c ... Accumulate the amounts flowed
      accum(UNLABL) = accum(UNLABL) + cflow1
      accum(LABELD) = accum(LABELD) + cflow2

      return
      end


text:
      subroutine sub1(num)
      write(6, *) 'sub1'
      if (num .eq. 1) then
        call sub2(num)
      endif
      end

      subroutine sub2(num)
      write(6, *) 'sub2'
      call sub1(num+1)
      end

      program prog
      write(6, *) 'start'
      call sub1(1)
      write(6, *) 'done'
      end


text:
module mod_monolis_direct_LU
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  implicit none

contains

!  subroutine  monolis_solver_direct_setup(monoPRM, monoCOM, monoMAT)
!    implicit none
!    type(monolis_prm) :: monoPRM
!    type(monolis_com) :: monoCOM
!    type(monolis_mat) :: monoMAT
!
!    call monolis_solver_direct_init(monoPRM, monoCOM, monoMAT)
!    call monolis_solver_direct_fact(monoPRM, monoCOM, monoMAT)
!  end subroutine monolis_solver_direct_setup
!
!  subroutine monolis_solver_direct_apply(monoPRM, monoCOM, monoMAT, X, Y)
!    implicit none
!    type(monolis_prm) :: monoPRM
!    type(monolis_com) :: monoCOM
!    type(monolis_mat) :: monoMAT
!    real(kind=kdouble) :: X(:), Y(:)
!
!    call monolis_solver_direct_solv(monoPRM, monoCOM, monoMAT)
!  end subroutine monolis_solver_direct_apply
!
!  subroutine monolis_solver_direct_clear(monoPRM, monoCOM, monoMAT)
!    implicit none
!    type(monolis_prm) :: monoPRM
!    type(monolis_com) :: monoCOM
!    type(monolis_mat) :: monoMAT
!
!    call monolis_solver_direct_clear(monoPRM, monoCOM, monoMAT)
!  end subroutine monolis_solver_direct_clear
!
!  subroutine monolis_solver_direct_init(monoPRM, monoCOM, monoMAT)
!    implicit none
!    type(monolis_prm) :: monoPRM
!    type(monolis_com) :: monoCOM
!    type(monolis_mat) :: monoMAT
!
!    if(monoCOM%commsize == 0) isEntire = .true.
!    call monolis_matrix_get_fillin(hecMESH, hecT, idxU, itemU, NPU)
!    call hecmw_matrix_copy_with_fillin(hecMESH, hecT, idxU, itemU, AU, NPU)
!  end subroutine monolis_solver_direct_init

end module mod_monolis_direct_LU

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
! KGEN-generated Fortran source file
!
! Filename    : rrlw_kg04.f90
! Generated at: 2015-07-06 23:28:45
! KGEN version: 0.4.13

    MODULE rrlw_kg04
        USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
        USE shr_kind_mod, ONLY: r8 => shr_kind_r8
        !      use parkind ,only : jpim, jprb
        IMPLICIT NONE
        !-----------------------------------------------------------------
        ! rrtmg_lw ORIGINAL abs. coefficients for interval 4
        ! band 4:  630-700 cm-1 (low - h2o,co2; high - o3,co2)
        !
        ! Initial version:  JJMorcrette, ECMWF, jul1998
        ! Revised: MJIacono, AER, jun2006
        !-----------------------------------------------------------------
        !
        !  name     type     purpose
        !  ----   : ----   : ---------------------------------------------
        !fracrefao: real
        !fracrefbo: real
        ! kao     : real
        ! kbo     : real
        ! selfrefo: real
        ! forrefo : real
        !-----------------------------------------------------------------
        !-----------------------------------------------------------------
        ! rrtmg_lw COMBINED abs. coefficients for interval 4
        ! band 4:  630-700 cm-1 (low - h2o,co2; high - o3,co2)
        !
        ! Initial version:  JJMorcrette, ECMWF, jul1998
        ! Revised: MJIacono, AER, jun2006
        !-----------------------------------------------------------------
        !
        !  name     type     purpose
        !  ----   : ----   : ---------------------------------------------
        ! absa    : real
        ! absb    : real
        !fracrefa : real
        !fracrefb : real
        ! ka      : real
        ! kb      : real
        ! selfref : real
        ! forref  : real
        !-----------------------------------------------------------------
        INTEGER, parameter :: ng4  = 14
        REAL(KIND=r8) :: fracrefa(ng4,9)
        REAL(KIND=r8) :: fracrefb(ng4,6)
        REAL(KIND=r8) :: absa(585,ng4)
        REAL(KIND=r8) :: absb(1175,ng4)
        REAL(KIND=r8) :: selfref(10,ng4)
        REAL(KIND=r8) :: forref(4,ng4)
        PUBLIC kgen_read_externs_rrlw_kg04
    CONTAINS

    ! write subroutines

    ! module extern variables

    SUBROUTINE kgen_read_externs_rrlw_kg04(kgen_unit)
        INTEGER, INTENT(IN) :: kgen_unit
        READ(UNIT=kgen_unit) fracrefa
        READ(UNIT=kgen_unit) fracrefb
        READ(UNIT=kgen_unit) absa
        READ(UNIT=kgen_unit) absb
        READ(UNIT=kgen_unit) selfref
        READ(UNIT=kgen_unit) forref
    END SUBROUTINE kgen_read_externs_rrlw_kg04

    END MODULE rrlw_kg04

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
subroutine evaluate_tidalfield(peps,ppos,ptide)
 include 'globals.h'
 real,intent(in) :: ppos(3),peps
 real,intent(inout) :: ptide(6)

 ptide(1:6)=0
 if(selfgrav) call system_tidalfield(peps,ppos,ptide)
 if(fixthalo) call halotidalfield(ppos,ptide)

end subroutine

subroutine system_tidalfield(peps,ppos,ptide)
 include 'globals.h'
 real,intent(in) :: ppos(3),peps
 real,intent(inout) :: ptide(6)
 integer i,nterms
! periodic -
  if(.NOT.directsum) then
   nterms=0
   call treewalk(root,ppos,peps,0.,0.,nterms)
  else
   nterms=nbodies
   do i=1,nbodies
    bodlist(i)=i
   enddo
  endif
  call tidalsum(peps,ppos,ptide,nterms)
end subroutine

subroutine tidalsum(peps,ppos,ptide,nterms)
 include 'globals.h'
 integer, intent(in):: nterms
 real,intent(in) :: peps,ppos(3)
 real,intent(inout):: ptide(6)
 real :: dx,dy,dz,dist2,dist,epseff,epseffi,drdeldrg, &
  &  rinveff,r3inveff,r5inveff,drsm,rhinv
 integer :: i,q,smindex
 do i=1,nterms
   q=bodlist(i)
   dx=ppos(1)-pos(q,1)
   dy=ppos(2)-pos(q,2)
   dz=ppos(3)-pos(q,3)
   dist2=dx*dx+dy*dy+dz*dz
   dist=SQRT(dist2)
   epseff=peps+epsgrav(q)
   if(dist.GE.epseff) then
    rinveff=1./dist
    r3inveff=rinveff**3
    r5inveff=rinveff**5
   else
    epseffi=1./epseff
    rhinv=dist*epseffi
    if(rhinv.LT.0.5) then
      r3inveff=4./3.+rhinv**2*(4.*rhinv-4.8)
      r5inveff=12*rhinv-9.6
    else
      r3inveff=8./3.-6.*rhinv+4.8*rhinv**2- &
            4./3.*rhinv**3-1./120./rhinv**3
      r5inveff=9.6-4.*rhinv-6./rhinv+1./40./rhinv**5
    endif
    r3inveff=r3inveff*8*epseffi**3
    r5inveff=-r5inveff*8./3.*epseffi**5
   endif

   ptide(1)=ptide(1)+mass(q)*(-r3inveff+3*r5inveff*dx**2)
   ptide(4)=ptide(4)+mass(q)*(-r3inveff+3*r5inveff*dy**2)
   ptide(6)=ptide(6)+mass(q)*(-r3inveff+3*r5inveff*dz**2)
   ptide(2)=ptide(2)+mass(q)*3*r5inveff*dx*dy
   ptide(3)=ptide(3)+mass(q)*3*r5inveff*dx*dz
   ptide(5)=ptide(5)+mass(q)*3*r5inveff*dy*dz

   if(usequad.AND.q.GE.nbods1) then
    call terror("tide+quad tbd")
   endif

 enddo
end subroutine tidalsum

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
! This subroutine calculates the absolute value of the determinant of a matrix
! using a LAPACK Q-R factorization. The absolute value of the determinant is
! the sum of the absolute values of the elements on the diagonal of the R
! matrix
! Calculate the absolute value of the determinant of a square NXN
! matrix using an LU factorization by LAPACK

subroutine abs_determinant(A, N, det)

   use params, only: dp

   integer N       !Size of the matrix A
   real(dp) A(N, N)
   real(dp) det
   integer LDA
   integer inFO
   integer I
   real(dp) LU(N, N)
   integer IPIV(N, N)

   LDA = N
   LWORK = N
   LU = A

   CALL DGETF2(N, N, A, LDA, IPIV, inFO)
   !Returned LU  matrix contains the U matrix from the LU factorization
   !on the upper diagonal. Take the absolute value of the product of the diagonals
   !to get the absolute value of the determinant

   det = 1.

   do I = 1, N
      det = det*A(I, I)
   ENDdo
   det = ABS(det)

END subroutine abs_determinant

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


text:
JeffSpeckles think Front Page this site has a lot of misinformation presented as facts because of the nature of the wiki:WikiPedia:Wiki medium. For instance, JeffSpeckles is an XastronautX XpaparazzoX XgigoloX penguin curator.  His advice:  Let the reader beware.

Note:  The user that posts under this account is not the Jeff Speckels that used to attend UCD:
http://www.sonicspike.net/pictures/albums/4/show/54

20070308 16:12:15 nbsp Jeff Speckels STRIKES AGAIN!!!!!! Users/JarrettNoble

20070319 12:20:11 nbsp If he doesnt use the wiki, he doesnt need a user page. Its like creating a facebook or myspace profile for a person who doesnt use either. Users/WilliamLewis

20070319 18:39:32 nbsp SaraHenry never created a user page and neither did DanteGarzon. You obviously dont understand the purpose of a userpage. Users/WilliamLewis

20070319 23:21:30 nbsp He doesnt qualify as a public figure and isnt notable for a public page.  Users/WilliamLewis

20070321 01:23:14 nbsp Hi Im William Lewis and since I keep deleting Dante and it keeps recurring, Im just gonna delete the comments from my site in order to act like there was never any conflict at all.  I delete, I however do not debate Users/JarrettNoble

20070703 20:23:22 nbsp robblerobble youre a jerk. robble. Users/ArlenAbraham

20071111 10:45:28 nbsp Was http://daviswiki.org/Fuzio?actiondiff&version268&version167 that not true? Users/EdWins

20110623 19:03:16 nbsp Hi Jeff Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! { dg-do compile }
! { dg-require-effective-target vect_double }
      subroutine shell(nx,ny,nz,q,dq)
      implicit none
      integer i,j,k,l,nx,ny,nz
      real*8 q(5,nx,ny),dq(5,nx,ny)
         do j=1,ny
            do i=1,nx
               do l=1,5
                  q(l,i,j)=q(l,i,j)+dq(l,i,j)
               enddo
            enddo
         enddo
      return
      end
! { dg-final { scan-tree-dump "vectorized 1 loops" "vect" { xfail { vect_no_align && { ! vect_hw_misalign } } } } }


text: Older Americans Month is a national event, but locally both the City of Davis City and the Senior Center http://www.cityofdavis.org/pcs/seniors/olderamericansmonth.cfm recognize it locally as well.


text:
      integer nproc
      common/nproc/nproc


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE VAC2
!                                                                      *
!   This routine sets up the second-order vacuum polarization poten-   *
!   tial using  equations (1) and (4) of  L Wayne Fullerton and  G A   *
!   Rinker, Jr,  Phys Rev A  13 (1976) 1283-1287.  The  potential is   *
!   accumulated  in  array  TB(I), I = 1, ..., N  which is in COMMON   *
!   block /TATB/ .                                                     *
!                                                                      *
!   Call(s) to: [LIB92]: QUAD.                                         *
!               [RCI92]: FUNK.                                         *
!                                                                      *
!   Written by Farid A Parpia, at Oxford    Last update: 09 Oct 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE def_C
      USE grid_C
      USE npar_C
      USE ncdist_C
      USE tatb_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE funk_I
      USE quad_I
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, K
      REAL(DOUBLE) :: EPSI, TWOCV, FACTOR, RI, X, TBI, RK, XK, XI, XM, XP
!-----------------------------------------------
!
!   Overall initialization
!
      EPSI = PRECIS*PRECIS
      TWOCV = CVAC + CVAC
!
!   Potential for a point nucleus: equation (1)
!   (this is also the asymptotoc form for a finite nucleus)
!
      FACTOR = -(2.0D00*Z)/(3.0D00*PI*CVAC)
!
      TB(1) = 0.0D00
!
      I = 1
    1 CONTINUE
      I = I + 1
!
      RI = R(I)
      X = TWOCV*RI
      TBI = (FACTOR/RI)*FUNK(X,1)
!
      IF (ABS(TBI) >= EPSI) THEN
         TB(I) = TBI
         IF (I < N) GO TO 1
      ELSE
         TB(I:N) = 0.0D00
      ENDIF
!
!   Potential for a finite nucleus: equation (4)
!
      IF (NPARM == 2) THEN
!
         FACTOR = -2.0D00/(3.0D00*CVAC**2)
!
!   Set up integrand
!
         TB(1) = 0.0D00
!
         K = 1
    3    CONTINUE
         K = K + 1
!
         RK = R(K)
         XK = TWOCV*RK
!
         TA(1) = 0.0D00
         DO I = 2, MTP
            XI = TWOCV*R(I)
            XM = ABS(XK - XI)
            XP = XK + XI
            TA(I) = (FUNK(XM,0) - FUNK(XP,0))*ZDIST(I)
         END DO
!
         CALL QUAD (X)
!
         X = X*FACTOR/RK
!
!   Get out of loop if the asymptotic value has been attained
!
         IF (ABS(X) >= EPSI) THEN
            IF (ABS((X - TB(K))/X) > 1.0D-05) THEN
               TB(K) = X
               IF (K < N) GO TO 3
            ENDIF
         ENDIF
!
      ENDIF
!
      RETURN
      END SUBROUTINE VAC2

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
      SUBROUTINE ALLMAG(MODEL,TM,RKM,ST,CT,SPH,CPH,BR,BT,BP,B)
C ****  GEOCENTRIC VERSION OF GEOMAGNETIC FIELD ROUTINE
C ****  SHORT DECK, USES SUBSCRIPTED VARIABLES AND DO LOOPS
C ****  EXECUTION TIME PER CALL 3 TIMES GREATER THAN LONG DECK
C ****  PROGRAM DESIGNED AND TESTED BY E G STASSINOPOULOS AND G D MEAD,
C ****  CODE 641, NASA GODDARD SPACE FLT CTR, GREENBELT, MD 20771
C  *****  INPUT  MODEL    CHOICE OF 7 MODELS - SEE BELOW
C  *****         RKM      GEOCENTRIC DISTANCE IN KILOMETERS
C  *****         TM       TIME IN YEARS FOR DESIRED FIELD
C  *****         ST,CT    SIN & COS OF GEOCENTRIC COLATITUDE
C  *****         SPH,CPH  SIN & COS OF EAST LONGITUDE
C  ***** OUTPUT  BR,BT,BP GEOCENTRIC FIELD COMPONENTS IN GAUSS
C  *****         B        FIELD MAGNITUDE IN GAUSS
C  *****  NOTE  FOR GREATEST EFFICIENCY, COMPLETE ALL CALCULATIONS WITH
C               ONE MODEL AND ONE TIME BEFORE CHANGING MODELS OR TIME.
      IMPLICIT REAL*8(A-H,O-Z)
C     REAL*8 LABEL(4,7) /  'HENDRICKS&CAIN 99-TERM GSFC 9/65 CAIN ET.AL.
C    A 120-TERM GSFC 12/66 CAIN&LANGEL 143-TERM POGO 10/68 CAIN&SWEENEY
C    B120-TERM POGO 8/69    IGRF 1965.0 80-TERM 10/68    LEATON MALIN EV
C    CANS 80-TERM 1965   HURWITZ US C&GS 168-TERM 1970'/
      DIMENSION T0(7),NMX(7),ISUM(7,3),G(13,13)
      DATA T0/4*1960.,2*1965.,1970./,NMX/10,11,12,11,9,9,13/
      INTEGER LSUM(7,3)/-1646106,-1795169,-1865298,-1777057,-158472,
     A-156856,-2191704,-62661,-96778,-181519,-83555,-9569,-9599,
     B-8593,1,-10618,5*1/
      INTEGER*4 G1(13,13),GT1(13,13),GTT1(13,13),G2(13,13),GT2(13,13),
     1 GTT2(13,13),G3(13,13),GT3(13,13),GTT3(13,13),G4(13,13),
     2 GT4(13,13),GTT4(13,13),G5(13,13),GT5(13,13),GTT5(13,13),
     3 G6(13,13),GT6(13,13),GTT6(13,13),G7(13,13),GT7(13,13),GTT7(13,13)
     4 ,LG(13,13,7),LGT(13,13,7),LGTT(13,13,7)
      REAL*4     GG(13,13,7),GGT(13,13,7),GGTT(13,13,7),SHMIT(13,13)
      EQUIVALENCE (G1(1),GG(1),LG(1)), (GT1(1),GGT(1),LGT(1)),
     A     (GTT1(1),GGTT(1),LGTT(1)),
     B (G2(1),LG(1,1,2)), (GT2(1),LGT(1,1,2)), (GTT2(1),LGTT(1,1,2)),
     C (G3(1),LG(1,1,3)), (GT3(1),LGT(1,1,3)), (GTT3(1),LGTT(1,1,3)),
     D (G4(1),LG(1,1,4)), (GT4(1),LGT(1,1,4)), (GTT4(1),LGTT(1,1,4)),
     E (G5(1),LG(1,1,5)), (GT5(1),LGT(1,1,5)), (GTT5(1),LGTT(1,1,5)),
     F (G6(1),LG(1,1,6)), (GT6(1),LGT(1,1,6)), (GTT6(1),LGTT(1,1,6)),
     G (G7(1),LG(1,1,7)), (GT7(1),LGT(1,1,7)), (GTT7(1),LGTT(1,1,7))
C  ***** THE FOLLOWING DATA CARDS CONTAIN THE FIELD COEFFICIENTS
C  ***** FOR THE FOLLOWING SEVEN MODELS
C  *****  G1,GT1  HENDRICKS & CAIN      99-TERM  GSFC  9/65  EPOCH 1960.
C  *****  G2,GT2,GTT2  CAIN ET. AL.    120-TERM  GSFC 12/66  EPOCH 1960.
C  *****  G3,GT3  CAIN & LANGEL        143-TERM  POGO 10/68  EPOCH 1960.
C  *****  G4,GT4  CAIN & SWEENEY       120-TERM  POGO  8/69  EPOCH 1960.
C  *****  G5,GT5  IGRF 1965.0           80-TERM       10/68  EPOCH 1965.
C  *****  G6,GT6  LEATON MALIN & EVANS 1965        80-TERM   EPOCH 1965.
C  *****          FOR MODEL 6 (LME 1965) SET RKM = 6371.2 + ALTITUDE
C  *****  G7,GT7  HURWITZ US COAST & GEODETIC S.  168-TERM   EPOCH 1970.
      DATA G1 / 10, -304249,-15361,13009,9576,-2277,498,709,48,99,3*0,
     A 57748,-21616,30002,-19870,8028,3595,607,-572,67,29,3*0,-19498,
     B 2043,15853,12904,5026,2313,45,56,-88,74,3*0,-4310,2308,-1300,8712
     C ,-3940,-312,-2417,75,-138,-156,3*0,1520,-2684,29,-2505,2714,
     D -1573,-12,-244,-33,114,3*0,86,1212,-1160,-1104,799,-652,5,-15,71,
     E 111,3*0,-119,1028,609,-272,-124,-116,-1091,141,-56,10,3*0,-540,
     F -244,-91,22,276,-211,-201,58,117,4*0,69,-122,58,-170,26,236,-25,
     G -160,64,16,3*0,-220,156,51,-35,-18,96,121,2,-25,15,42*0 /
      DATA GT1 / 100, 2059,-2907,266,-86,255,-70,6*0,-394,602,121,-1003,
     H 194,-8,99,6*0,-1369,-1578,-70,163,-117,153,85,6*0,649,293,-924,
     I -130,-54,-42,211,6*0,-177,-154,318,-548,-417,-72,157,6*0,304,288,
     J -186,125,80,164,-9,6*0,-139,12,153,-73,-6,45,6,84*0/
      DATA GTT1 /1,168*0/
      DATA G2 / 10, -304012,-15401,13071,9493,-2335,492,722,85,104,-29,
     A 2*0,57782,-21638,29979,-19889,8035,3557,575,-537,65,58,-9,
     B 2*0,-19320,2029,15903,12768,5029,2284,-8,79,-93,75,-22,2*0,-4254,
     C 2278,-1338,8812,-3977,-288,-2383,156,-96,-151,8,2*0,1603,-2743,
     D 23,-2466,2665,-1579,-15,-243,-61,121,-28,2*0,51,1178,-1148,-1089,
     E 824,-622,-20,-36,55,47,64,2*0,-121,1044,566,-234,-148,-133,-1089,
     F 155,-81,2,47,2*0,-537,-274,-81,70,243,-225,-214,36,130,16,-2,2*0,
     G 54,-117,42,-153,46,219,-7,-171,74,9,18,2*0,-224,138,63,-30,-19,
     H 90,115,1,-15,2,20,2*0,-1,45,-10,26,-44,-13,-36,40,10,-20,11,28*0/
      DATA GT2 / 100, 1403,-2329,-93,145,161,-42,-57,35,-10,-1,2*0,-371,
     I 876,-9,-1062,90,60,82,-34,50,-13,-13,2*0,-1431,-1662,-456,231,
     J -175,334,82,-144,170,-120,88,2*0,520,253,-698,-589,66,-4,235,-90,
     K -11,8,-18,2*0,-219,-14,188,-652,-301,-60,83,3,34,-8,17,2*0,224,
     L 159,-261,50,-12,176,1,-60,-7,-39,-2,2*0,5,9,255,-119,33,84,23,-17
     M ,43,-36,5,2*0,-96,1,43,75,-33,49,90,-64,-15,47,17,2*0,-50,-21,3,
     N -79,5,10,-36,-43,-42,37,16,2*0,66,54,3,35,-3,-1,45,-5,75,-46,31,
     O 2*0,-61,-64,2,5,-63,-7,7,-3,-2,-45,-23,28*0/
      DATA GTT2 /1000,-62,-154,-123,1,45,-6,-14,6,-5,-3,2*0,-43,114,-18,
     P-27,-44,1,15,-6,8,-1,-3,2*0,54,-16,-253,28,17,75,10,-34,39,-27,20,
     Q 2*0,95,-7,79,-183,7,8,50,-4,-8,5,-8,2*0,4,56,-35,-47,-97,15,-11,
     R -6,15,-7,7,2*0,-46,7,-7,1,-24,56,26,-27,-2,-6,1,2*0,20,-11,15,
     S -29,29,-10,23,-1,5,-9,1,2*0,-14,16,14,5,-8,16,11,-4,-8,6,1,2*0,
     T -15,-12,5,-11,0,-3,-9,-3,-7,5,5,2*0,22,7,-2,9,6,-1,9,-4,19,-9,4,
     U 2*0,-12,-14,1,1,-11,-1,1,-1,1,-6,-2,28*0/
      DATA G3 / 10, -304650,-15414,13258,9591,-2343,491,759,74,110,-26,
     A 23,0,57910,-21633,29763,-19837,8196,3577,545,-524,60,66,-20,-18,
     B 0,-19772,1566,16075,13169,4864,2339,48,80,-81,18,10,-21,0,-4453,
     C 2334,-949,8420,-3724,-210,-2491,100,-92,-125,-55,55,0,1354,-2667,
     D 207,-2415,2562,-1471,17,-367,-8,158,-7,-15,0,169,1133,-1287,-1151
     E ,1303,-452,-37,-83,91,17,75,24,0,-96,1064,568,-272,-149,-43,-916,
     F 66,-114,26,78,-35,0,-579,-250,-8,63,95,-117,-376,-227,79,87,17,
     G -13,0,101,-130,115,-164,55,223,-49,-262,351,51,-53,25,0,-204,144,
     H 6,-15,14,34,148,24,-9,-24,13,-12,0,11,9,-3,75,-23,14,-5,43,80,
     I -137,-27,127,0,-8,44,-1,-39,-6,18,-32,8,-59,-17,105,50,14*0/
      DATA GT3 / 100,2542,-2390,-559,-62,272,-61,-89,61,-24,-1,3,0,-466,
     J 988,350,-1152,-251,48,106,-21,-12,30,-9,11,0,-707,-1070,-214,-441
     K ,-122,317,62,-108,87,4,12,5,0,848,68,-1489,287,-296,-246,396,70,
     L -33,4,19,-30,0,345,-39,-87,-652,86,-89,-94,107,-14,-40,-20,1,0,5,
     M 300,32,311,-635,-315,149,96,-85,-28,-2,-34,0,-26,-48,258,-80,50,
     N 82,-167,101,99,-57,-43,48,0,-87,-46,-102,25,188,-243,232,523,81,
     O -132,-33,52,0,-15,-10,-122,-26,15,-37,29,91,-498,-14,103,-19,0,
     P -38,16,67,-14,-83,130,-33,-38,99,50,22,-3,0,21,5,54,-26,-30,-3,
     Q -39,-2,-104,79,46,-165,0,35,-26,-17,17,18,-50,23,-34,37,22,-155,
     R -40,14*0/
      DATA GTT3 /1,168*0/
      DATA G4 / 10,-304708,-15425,13334,9647,-2375,448,793,99,96,-17,
     A 2*0,57571,-21702,29893,-19826,8108,3566,594,-516,32,93,-22,2*0,
     B -19793,2661,15559,12922,5068,2498,-37,-3,-56,31,13,2*0,-4249,
     C 2417,-1740,8336,-3978,-143,-2324,89,-165,-120,16,2*0,1344,-3037,
     D 194,-2764,2247,-1497,96,-335,-33,153,-22,2*0,51,1080,-1073,-1083,
     E 1171,-757,20,-33,50,7,94,2*0,-76,1181,583,-181,-270,1,-831,100,
     F -120,8,87,2*0,-544,-212,-87,55,151,-236,-278,39,102,4,3,2*0,98,
     G -162,99,-189,106,206,-2,-207,187,62,-24,2*0,-254,128,31,-25,-21,
     H 73,127,47,7,-38,-1,2*0,29,35,-7,66,-50,10,-28,21,42,-88,53,28*0/
      DATA GT4 / 100,2682,-2366,-724,-157,359,12,-160,19,17,-3,2*0,225,
     I 1003,150,-1142,-118,58,38,-26,27,-8,-8,2*0,-684,-2832,792,84,
     J -536,-27,235,72,33,-46,17,2*0,449,-96,177,327,102,-326,128,86,83,
     K -9,-87,2*0,369,564,-109,-205,834,-108,-277,84,42,-37,-12,2*0,234,
     L 401,-424,63,-503,504,8,-57,0,-3,-33,2*0,-65,-238,249,-170,234,
     M -259,-130,101,49,-48,-33,2*0,-168,-114,58,123,94,40,60,-140,73,
     N 54,-21,2*0,1,39,-106,-9,-49,56,-67,-8,-148,-13,27,2*0,48,42,17,
     O -41,-22,21,1,-113,16,33,49,2*0,-14,-37,51,-2,4,-19,7,40,-53,31,
     P -75,28*0/
      DATA GTT4 /1,168*0/
      DATA G5 / 1, -30339,-1654,1297,958,-223,47,71,10,4*0,5758,-2123,
     A 2994,-2036,805,357,60,-54,9,4*0,-2006,130,1567,1289,492,246,4,0,
     B -3,4*0,-403,242,-176,843,-392,-26,-229,12,-12,4*0,149,-280,8,-265
     C ,256,-161,3,-25,-4,4*0,16,125,-123,-107,77,-51,-4,-9,7,4*0,-14,
     D 106,68,-32,-10,-13,-112,13,-5,4*0,-57,-27,-8,9,23,-19,-17,-2,12,
     E 4*0,3,-13,5,-17,4,22,-3,-16,6,56*0/
      DATA GT5 / 10, 153,-244,2,-7,19,-1,-5,1,4*0,-23,87,3,-108,2,11,-3,
     F -3,4,4*0,-118,-167,-16,7,-30,29,11,-7,6,4*0,42,7,-77,-38,-1,6,19,
     G -5,5*0,-1,16,29,-42,-21,0,-4,3,5*0,23,17,-24,8,-3,13,-4,0,-1,4*0,
     H-9,-4,20,-11,1,9,-2,-2,3,4*0,-11,3,4,2,4,2,3,-6,-3,4*0,1,-2,-3,-2,
     I -3,-4,-3,-3,-5,56*0/
      DATA GTT5 /1,168*0/
      DATA G6 / 1, -30375,-1648,1164,930,-179,42,77,11,4*0,5769,-2087,
     A 2954,-2033,811,357,55,-56,23,4*0,-1995,116,1579,1299,490,248,12,
     B 8,-6,4*0,-389,230,-141,880,-402,-20,-239,5,-17,4*0,142,-276,5,
     C -264,262,-171,16,-35,5,4*0,30,135,-123,-100,84,-64,8,-16,20,4*0,
     D -18,101,60,-32,-27,-12,-110,9,-1,4*0,-47,-35,-9,2,27,-17,-24,2,
     E 12,4*0,5,-7,3,-20,8,26,10,-12,7,56*0/
      DATA GT6 / 10, 155,-266,0,6,8,7*0,6,83,-13,-95,10,4,-5,6*0,-114,
     F -182,13,-19,-22,16,18,6*0,32,16,-85,-6,2,-3,14,6*0,30,-7,27,-27,
     G -30,-11,6,6*0,19,23,-18,14,5,17,2,6*0,-22,2,9,-21,-1,-2,-22,84*0/
      DATA GTT6 /1,168*0/
      DATA G7/10,-302059,-17917,12899,9475,-2145,460,734,121,107,-39,16,
     A -4,57446,-20664,29971,-20708,8009,3595,651,-546,77,57,-26,-31,30,
     B -20582,430,16086,12760,4579,2490,95,46,-32,23,7,-36,5,-3699,2456,
     C -1880,8334,-3960,-290,-2188,175,-124,-110,-19,37,-3,1617,-2758,
     D 185,-2788,2436,-1669,20,-210,-44,131,-15,-3,-13,157,1420,-1310,
     E -911,808,-582,-22,-32,45,33,74,-6,4,-171,1146,625,-323,-78,38,
     F -1125,143,34,2,46,-8,-14,-666,-265,-34,81,209,-240,-186,41,125,
     G 15,6,1,-12,121,-160,22,-176,46,189,-46,-187,94,9,-8,2,-12,-174,
     H 163,14,-27,-32,80,137,-4,-14,-4,22,-24,-1,27,19,0,35,-45,22,-31,
     I 56,-1,-63,14,4,10,-2,26,-26,-9,21,-1,18,-14,-28,-17,-14,6,-4,-3,
     J 4,9,-1,-10,26,-32,13,-6,-19,7,19,12/
      DATA GT7/10,231,-244,-19,-7,12,-7,0,3,4*0,-46,112,-1,-90,-6,7,6,
     K -3,3,4*0,-104,-166,40,-20,-36,12,14,3,4,4*0,72,21,-52,-54,-11,0,
     L 17,6,1,4*0,22,-5,14,-24,-23,-15,6,3,-1,4*0,1,25,-14,9,1,11,-3,2,
     M -3,4*0,-5,11,2,-3,7,22,-5,1,9,4*0,-17,-3,7,1,-2,-3,-2,-1,-2,4*0,
     N 2,-6,-3,-4,1,-2,-2,-1,6,56*0/
      DATA GTT7 /1,168*0/
      DATA SHMIT(1,1) / 0.0 /, TMOLD / 0.0 /, MODOLD / 0 /
C  *****  SUBSCRIPTED DO-LOOP VERSION BEGINS HERE
      DIMENSION CONST(13,13),FN(13),FM(13)
      DIMENSION P(13,13),DP(13,13),SP(13),CP(13)
      DATA P(1,1),CP(1),DP(1,1),SP(1) / 2*1.,2*0. /
C  *****  BEGIN PROGRAM
      IF(SHMIT(1,1).EQ.-1.)   GO TO 8
C  *****  INITIALIZE * ONCE ONLY, FIRST TIME SUBROUTINE IS CALLED
      SHMIT(1,1)=-1.
      DO 18 N=1,13
      FN(N)=N
      DO 18 M=1,13
      FM(M)=M-1
   18 CONST(N,M) = FLOAT((N-2)**2-(M-1)**2) / ((2*N-3)*(2*N-5))
      DO 2 N=2,13
      SHMIT(N,1) = (2*N-3) * SHMIT(N-1,1) / (N-1)
      JJ=2
      DO 2 M=2,N
      SHMIT(N,M) = SHMIT(N,M-1) * SQRT(FLOAT((N-M+1)*JJ)/(N+M-2))
      SHMIT(M-1,N)=SHMIT(N,M)
    2 JJ = 1
      DO 7 K=1,7
      F1=LG(1,1,K)
      F2=LGT(1,1,K)
      F3=LGTT(1,1,K)
      NMAX=NMX(K)
      L = 0
      DO 3 I=1,3
    3 ISUM(K,I) = 0
      DO 4 N=1,NMAX
      DO 4 M=1,NMAX
      L = L+1
      ISUM(K,1)=ISUM(K,1)+L*LG(N,M,K)
      ISUM(K,2)=ISUM(K,2)+L*LGT(N,M,K)
    4 ISUM(K,3)=ISUM(K,3)+L*LGTT(N,M,K)
      DO 6 I=1,3
      IF(ISUM(K,I).EQ.LSUM(K,I))  GO TO 6
C  *****  ERROR IN DATA CARDS - NOTE WRITE AND STOP STATEMENTS
      PRINT 5,   K,I,LSUM(K,I),ISUM(K,I)
    5 FORMAT(///29H DATA WRONG IN ALLMAG--MODEL ,I2,3X,2HI=,I1,3X,
     A17HPRECALCULATED SUM,I10,3X,17HTHIS MACHINE GETS,I10)
      STOP
    6 CONTINUE
      DO 7 N=1,NMAX
      DO 7 M=1,NMAX
      GG(N,M,K)=LG(N,M,K)*SHMIT(N,M)/F1
      GGT(N,M,K)=LGT(N,M,K)*SHMIT(N,M)/F2
    7 GGTT(N,M,K)=LGTT(N,M,K)*SHMIT(N,M)/F3
    8 IF((MODEL.EQ.MODOLD).AND.(TM.EQ.TMOLD))  GO TO 11
C  *****  NOTE WRITE STATEMENT - NEW MODEL OR NEW TIME
C    DO NOT WRITE
C     PRINT 9,   MODEL,(LABEL(I,MODEL),I=1,4),TM
C   9 FORMAT('0 MODEL USED IS NUMBER',I2,2X,4A8,'  FOR TM =',F9.3/)
      IF(MODEL.LT.1.OR.MODEL.GT.7) STOP
      MODOLD=MODEL
      TMOLD=TM
      NMAX=NMX(MODEL)
      T=TM-T0(MODEL)
      DO 10 N=1,NMAX
      DO 10 M=1,NMAX
   10 G(N,M)=GG(N,M,MODEL)+T*(GGT(N,M,MODEL)+GGTT(N,M,MODEL)*T)
C  *****  CALCULATION USUALLY BEGINS HERE
   11 SP(2)=SPH
      CP(2)=CPH
      DO 12 M=3,NMAX
      SP(M)=SP(2)*CP(M-1)+CP(2)*SP(M-1)
   12 CP(M)=CP(2)*CP(M-1)-SP(2)*SP(M-1)
      AOR=6371.2/RKM
      AR=AOR**2
      BR=0.0
      BT=0.0
      BP=0.0
      DO 17 N=2,NMAX
      AR=AOR*AR
      DO 17 M=1,N
      IF(M.EQ.N) GO TO 13
      P(N,M)=CT*P(N-1,M)-CONST(N,M)*P(N-2,M)
      DP(N,M)=CT*DP(N-1,M)-ST*P(N-1,M)-CONST(N,M)*DP(N-2,M)
      GO TO 14
   13 P(N,N)=ST*P(N-1,N-1)
      DP(N,N)=ST*DP(N-1,N-1)+CT*P(N-1,N-1)
   14 PAR=P(N,M)*AR
      IF(M.EQ.1) GO TO 15
      TEMP=G(N,M)*CP(M)+G(M-1,N)*SP(M)
      BP=BP-(G(N,M)*SP(M)-G(M-1,N)*CP(M))*FM(M)*PAR
      GO TO 16
   15 TEMP = G(N,M)
   16 BR=BR-TEMP*FN(N)*PAR
   17 BT=BT+TEMP*DP(N,M)*AR
    1 BR = BR / 100000.
      BT = BT / 100000.
      BP = BP / ST / 100000.
      B = DSQRT(BR*BR+BT*BT+BP*BP )
      RETURN
      END

ext: for


text:
module cfdio

  use mpi

  implicit none

contains

subroutine writedatafiles(psi, m, n, scale, comm)

  integer :: m, n, scale, comm
  double precision ::  psi(0:m+1, 0:n+1)

  double precision, allocatable :: vel(:,:,:)
  integer, allocatable :: rgb(:,:,:)

  double precision :: modvsq, hue
  integer :: size, rank, irank, i, j, k,  ierr, ix, iy

  integer, parameter :: tag = 1, iounitvel = 10, iounitcol = 11

  integer, dimension(MPI_STATUS_SIZE) :: status

  call mpi_comm_size(comm, size, ierr)
  call mpi_comm_rank(comm, rank, ierr)

! Compute local velocities and colours

  allocate(rgb(3,m,n))
  allocate(vel(2,m,n))

  do i = 1, m
     do j = 1, n

        vel(1,i,j) =   (psi(i,j+1)-psi(i,j-1)) / 2.0
        vel(2,i,j) = - (psi(i+1,j)-psi(i-1,j)) / 2.0

        modvsq = vel(1,i,j)**2 + vel(2,i,j)**2
        hue = modvsq**0.4

        call hue2rgb(hue, rgb(1,i,j), rgb(2,i,j), rgb(3,i,j))

     end do
  end do

!  Receive data

  if (rank == 0) then

     open(unit=iounitcol, file='colourmap.dat', form='formatted')
     open(unit=iounitvel, file='velocity.dat',  form='formatted')

     do irank = 0, size-1

        if (irank /= 0) then

           call mpi_recv(rgb, 3*m*n, MPI_INTEGER, irank, tag, &
                         comm, status, ierr)

           call mpi_recv(vel, 2*m*n, MPI_DOUBLE_PRECISION, irank, tag, &
                         comm, status, ierr)

        end if

!  Write out

        do j = 1, n

           iy = irank*n + j

           do i = 1, m

              ix = i

!  Write colour map of velocity magnitude at every point

              write(iounitcol,fmt='(i4,1x,i4,1x,i3,1x,i3,1x,i3)') &
                    ix, iy, rgb(1,i,j), rgb(2,i,j), rgb(3,i,j)

!  Only write velocity vectors every "scale" points

              if (mod(ix-1,scale) == (scale-1)/2 .and. &
                  mod(iy-1,scale) == (scale-1)/2         ) then

                 write(iounitvel,fmt='(i4,1x,i4,1x,g12.5,1x,g12.5)') &
                       ix, iy, vel(1,i,j), vel(2,i,j)
              end if

           end do
        end do
     end do

     close(unit=iounitcol)
     close(unit=iounitvel)

  else

     call mpi_ssend(rgb, 3*m*n, MPI_INTEGER, 0, tag, comm, ierr)
     call mpi_ssend(vel, 2*m*n, MPI_DOUBLE_PRECISION, 0, tag, comm, ierr)

  end if

end subroutine writedatafiles

subroutine writeplotfile(m, n, scale)

  integer :: m, n, scale
  integer, parameter :: iounit = 10

  open(unit=iounit, file='cfd.plt', form='formatted')

  write(iounit,*) 'set size square'
  write(iounit,*) 'set key off'
  write(iounit,*) 'unset xtics'
  write(iounit,*) 'unset ytics'

  write(iounit,fmt='('' set xrange ['',i4,'':'',i4, '']'')') 1-scale, m+scale
  write(iounit,fmt='('' set yrange ['',i4,'':'',i4, '']'')') 1-scale, n+scale

  write(iounit,fmt='('' plot "colourmap.dat" w rgbimage, "velocity.dat" u 1:2:&
       &('',i2,''*0.75*$3/sqrt($3**2+$4**2)):&
       &('',i2,''*0.75*$4/sqrt($3**2+$4**2)) &
       &with vectors  lc rgb "#7F7F7F"'')') scale, scale

  close(unit=iounit)

end subroutine writeplotfile

subroutine hue2rgb(hue, r, g, b)

  double precision :: hue

  integer :: r, g, b
  integer, parameter :: rgbmax = 255

  r = rgbmax*colfunc(hue-1.0)
  g = rgbmax*colfunc(hue-0.5)
  b = rgbmax*colfunc(hue    )

end subroutine hue2rgb

double precision function colfunc(x)

  double precision :: x, absx, val

  double precision, parameter :: x1 = 0.2, x2 = 0.5

  absx = abs(x)

  if (absx .gt. x2) then
     val = 0.0
  else if (absx .lt. x1) then
     val = 1.0
  else
     val = 1.0 - ((absx-x1)/(x2-x1))**2
  end if

  colfunc = val

end function colfunc

end module cfdio

max_stars_repo_name: davidhenty/CFD

max_issues_repo_name: davidhenty/CFD

max_forks_repo_name: davidhenty/CFD


text:
MODULE parameter_mod

  IMPLICIT NONE
  SAVE

  INTEGER, PARAMETER :: dp = selected_real_kind(14,300)

  INTEGER, PARAMETER :: nx=4096, ny=2048, nz=1
  INTEGER, PARAMETER :: nx1=nx-1, ny1=ny-1, nz1=nz-1
  INTEGER, PARAMETER :: nx2=nx-2
  INTEGER, PARAMETER :: nx3=nx-3
  INTEGER, PARAMETER :: nx4=nx-4
  INTEGER, PARAMETER :: n2x=nx/2, n2y=ny/2
! n2z=nz/2 computed in pstartup (so to switch 2D/3D cases)

  INTEGER, PARAMETER :: l_rst=500, l_cpu=501, l_div=502, l_Eng=503 , &
                        l_U=504, l_DPJ=505, l_EB=506, l_Gi=507

  INTEGER, PARAMETER :: nprocx = 20000
  INTEGER :: l_EB_me, l_U_me, l_DPJ_me, l_Gi_me, l_rst_me, n2z
  REAL(dp)    :: pgreco, ninvy, ninvz

  CHARACTER*4 :: prefix

  CHARACTER(LEN=2), PARAMETER :: space_dim = '2D'

END MODULE parameter_mod

MODULE box_mod

  use parameter_mod, only: dp, nx, ny, nz

  IMPLICIT NONE
  SAVE

  integer :: istart, ioutt, ioutx, nwx, nwy, nwz, ibc
  integer :: rk_ord, nx_lmbd, hall_on, divPe_on, flr_on

  REAL(dp)  :: rapm, de, de2, de12inv, gam, gm1, eta, ampl, alpha
  REAL(dp)  :: dt, ab1, ab2, ab3
  REAL(dp)  :: xl, yl, zl, dx, dy, dz, dxyz, hyl, hzl
  REAL(dp)  :: Tot_inv_ny, Tot_inv_nz, Tot_inv_nxnynz
  REAL(dp)  :: tempo, tmax, tt_last, tt_w, tx_last, tx_w
  REAL(dp)  :: nx_lambda, terzo, grdlmbd

  REAL(dp),DIMENSION (nx) :: x
  REAL(dp),DIMENSION (ny) :: y
  REAL(dp),DIMENSION (nz) :: z
  REAL(dp),DIMENSION (nx) :: lambda

END MODULE box_mod

MODULE fields_UJ_mod

  use parameter_mod, only: dp, nx, ny, nz

  IMPLICIT NONE
  SAVE

  REAL(dp) :: fac1_i, fac2_i, fac3_i
  REAL(dp) :: fac1_e, fac2_e, fac3_e

  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Uex, Uey, Uez, Uix, Uiy, Uiz
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Ux, Uy, Uz, Jx, Jy, Jz
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Tracciante
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Ti_para, Ti_perp, Te_para, Te_perp

CONTAINS

  SUBROUTINE allocate_fields_UJ( nxl, nyl, nzl )
     INTEGER :: nxl, nyl, nzl
     ALLOCATE( Uex( nxl, nyl, nzl ) )
     ALLOCATE( Uey( nxl, nyl, nzl ) )
     ALLOCATE( Uez( nxl, nyl, nzl ) )
     ALLOCATE( Uix( nxl, nyl, nzl ) )
     ALLOCATE( Uiy( nxl, nyl, nzl ) )
     ALLOCATE( Uiz( nxl, nyl, nzl ) )
     ALLOCATE( Ux( nxl, nyl, nzl ) )
     ALLOCATE( Uy( nxl, nyl, nzl ) )
     ALLOCATE( Uz( nxl, nyl, nzl ) )
     ALLOCATE( Jx( nxl, nyl, nzl ) )
     ALLOCATE( Jy( nxl, nyl, nzl ) )
     ALLOCATE( Jz( nxl, nyl, nzl ) )
     ALLOCATE( Tracciante( nxl, nyl, nzl ) )
     ALLOCATE( Ti_para( nxl, nyl, nzl ) )
     ALLOCATE( Ti_perp( nxl, nyl, nzl ) )
     ALLOCATE( Te_para( nxl, nyl, nzl ) )
     ALLOCATE( Te_perp( nxl, nyl, nzl ) )
     RETURN
  END SUBROUTINE

END MODULE fields_UJ_mod

MODULE fields_DP_mod

  use parameter_mod, only: dp, nx, ny, nz

  IMPLICIT NONE
  SAVE

  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: ax, ay, az, a_inv, c_s, c_s_inv
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: s, f, f_inv
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: alpha_1  , alpha_2  , segno
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: L_s_m, L_a_m, L_f_m
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: L_s_p, L_a_p, L_f_p
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: L_0
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Den, Dinv
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: pe_para, pe_perp, pi_para, pi_perp
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: nU_x, nU_y, nU_z
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: qi_para, qi_perp, qe_para, qe_perp
  REAL(dp),ALLOCATABLE, DIMENSION (:,:,:) :: Gi_xx, Gi_xy, Gi_xz, Gi_yz

CONTAINS

  SUBROUTINE allocate_fields_DP( nxl, nyl, nzl )
     INTEGER :: nxl, nyl, nzl
     ALLOCATE( ax( 2, nyl, nzl ) )
     ALLOCATE( ay( 2, nyl, nzl ) )
     ALLOCATE( az( 2, nyl, nzl ) )
     ALLOCATE( a_inv( 2, nyl, nzl ) )
     ALLOCATE( c_s( 2, nyl, nzl ) )
     ALLOCATE( c_s_inv( 2, nyl, nzl ) )
     ALLOCATE( s( 2, nyl, nzl ) )
     ALLOCATE( f( 2, nyl, nzl ) )
     ALLOCATE( f_inv( 2, nyl, nzl ) )
     ALLOCATE( alpha_1( 2, nyl, nzl ) )
     ALLOCATE( alpha_2( 2, nyl, nzl ) )
     ALLOCATE( segno( 2, nyl, nzl ) )
     ALLOCATE( L_s_m( 2, nyl, nzl ) )
     ALLOCATE( L_a_m( 2, nyl, nzl ) )
     ALLOCATE( L_f_m( 2, nyl, nzl ) )
     ALLOCATE( L_s_p( 2, nyl, nzl ) )
     ALLOCATE( L_a_p( 2, nyl, nzl ) )
     ALLOCATE( L_f_p( 2, nyl, nzl ) )
     ALLOCATE( L_0( 2, nyl, nzl ) )
     ALLOCATE( Den( nxl, nyl, nzl ) )
     ALLOCATE( Dinv( nxl, nyl, nzl ) )
     ALLOCATE( pi_para( nxl, nyl, nzl ) )
     ALLOCATE( pi_perp( nxl, nyl, nzl ) )
     ALLOCATE( pe_para( nxl, nyl, nzl ) )
     ALLOCATE( pe_perp( nxl, nyl, nzl ) )
     ALLOCATE( nU_x( nxl, nyl, nzl ) )
     ALLOCATE( nU_y( nxl, nyl, nzl ) )
     ALLOCATE( nU_z( nxl, nyl, nzl ) )
     ALLOCATE( qi_para( nxl, nyl, nzl ) )
     ALLOCATE( qi_perp( nxl, nyl, nzl ) )
     ALLOCATE( qe_para( nxl, nyl, nzl ) )
     ALLOCATE( qe_perp( nxl, nyl, nzl ) )
     ALLOCATE( Gi_xx( nxl, nyl, nzl ) )
     ALLOCATE( Gi_xy( nxl, nyl, nzl ) )
     ALLOCATE( Gi_xz( nxl, nyl, nzl ) )
     ALLOCATE( Gi_yz( nxl, nyl, nzl ) )
     RETURN
  END SUBROUTINE

END MODULE fields_DP_mod

MODULE fields_EB_mod

  use parameter_mod, only: dp, nx, ny, nz

  IMPLICIT NONE
  SAVE

   !REAL(dp) :: Bx0, By0, Bz0
   REAL(dp) :: B00, angolo, Aeq, Leq, DeltaDen

  REAL(dp), ALLOCATABLE, DIMENSION (:,:,:) :: Ex, Ey, Ez, Bx, By, Bz

CONTAINS

  SUBROUTINE allocate_fields_EB( nxl, nyl, nzl )
    INTEGER :: nxl, nyl, nzl
    ALLOCATE( Ex( nxl, nyl, nzl ) )
    ALLOCATE( Ey( nxl, nyl, nzl ) )
    ALLOCATE( Ez( nxl, nyl, nzl ) )
    ALLOCATE( Bx( nxl, nyl, nzl ) )
    ALLOCATE( By( nxl, nyl, nzl ) )
    ALLOCATE( Bz( nxl, nyl, nzl ) )
    RETURN
  END SUBROUTINE

END MODULE fields_EB_mod

MODULE deriv_mod

  use parameter_mod, only: dp, nx, ny, nz, nx1, nx2, n2z

  IMPLICIT NONE
  SAVE

! Derivate open (lungo x)

CHARACTER*1 TRANS

integer :: nrhs, ndb, info

integer :: ipv_d(nx), ipv_w(nx), ipv_q(nx)

real(dp) :: a_1_3x3, a_1_3x5, b_1_3x5
real(dp) :: a_1_bc, b_1_bc, c_1_bc, d_1_bc

REAL(dp),DIMENSION (nx)  :: dd_1
REAL(dp),DIMENSION (nx1) :: dm_1, dp_1
REAL(dp),DIMENSION (nx2) :: dw_1

REAL(dp),DIMENSION (nx)  :: wd_1
REAL(dp),DIMENSION (nx1) :: wm_1, wp_1
REAL(dp),DIMENSION (nx2) :: ww_1

REAL(dp),DIMENSION (nx)  :: qd_1
REAL(dp),DIMENSION (nx1) :: qm_1, qp_1
REAL(dp),DIMENSION (nx2) :: qq_1

! Derivate periodiche (lungo y,z)

REAL(dp):: ky_1(ny/2 - 1), ky_4(ny/2 - 1), savey(2 * ny + 15)
REAL(dp):: savez(2 * nz + 15)
REAL(dp), ALLOCATABLE, DIMENSION (:) :: kz_1, kz_4

END MODULE deriv_mod

MODULE filtro_per_mod

  use parameter_mod, only: dp, nx, ny, nz

  IMPLICIT NONE
  SAVE

real(dp) :: omega_x, omega_y, omega_z

!perche'? omega_x def anche qui????? refuso double-periodic??

REAL(dp):: work4x(nx/2), work4y(ny/2)
REAL(dp), ALLOCATABLE, DIMENSION (:) :: work4z

END MODULE filtro_per_mod

MODULE filtro_open_mod

  use parameter_mod, only: dp, nx, nx1, nx2

  IMPLICIT NONE
  SAVE

CHARACTER*1 TARNS

integer :: lrhs, ldb, info

integer :: ipv_f(nx)

real(dp) :: omega_x
real(dp) :: a_filt, b_filt, c_filt, d_filt
real(dp) :: a11_filt, a12_filt, a13_filt, a14_filt, a15_filt
real(dp) :: a21_filt, a22_filt, a23_filt, a24_filt, a25_filt
real(dp) :: a31_filt, a32_filt, a33_filt, a34_filt, a35_filt

REAL(dp),DIMENSION (nx)  :: filt_d
REAL(dp),DIMENSION (nx1) :: filt_m, filt_p
REAL(dp),DIMENSION (nx2) :: filt_2

END MODULE filtro_open_mod

MODULE poisson_open
  use parameter_mod, only: dp, nx, ny, nz, nx1, nx2, n2y, n2z

  IMPLICIT NONE
  SAVE
CHARACTER*1 TRNAS
integer :: irhs, idb, ifno
REAL(dp),DIMENSION (nx2) :: poiss_w

  integer , ALLOCATABLE, DIMENSION (:,:,:) :: ipv_p
  REAL(dp), ALLOCATABLE, DIMENSION (:,:)   :: pss_rhs_0, pss_rhs_1
  REAL(dp), ALLOCATABLE, DIMENSION (:,:,:) :: poiss_d
  REAL(dp), ALLOCATABLE, DIMENSION (:,:,:) :: poiss_m, poiss_p

CONTAINS

  SUBROUTINE allocate_poisson( n2z )
    INTEGER :: n2z

    ALLOCATE( ipv_p(nx,0:n2y-1,0:n2z-1) )
    ALLOCATE( pss_rhs_0(0:n2y-1,0:n2z-1) )
    ALLOCATE( pss_rhs_1(0:n2y-1,0:n2z-1) )
    ALLOCATE( poiss_d(nx,0:n2y-1,0:n2z-1) )
    ALLOCATE( poiss_m(nx-1,0:n2y-1,0:n2z-1) )
    ALLOCATE( poiss_p(nx-1,0:n2y-1,0:n2z-1) )

    RETURN
  END SUBROUTINE

END MODULE poisson_open

MODULE dom_distr_mod

  use parameter_mod, only: nprocx

  IMPLICIT NONE
  SAVE

  INTEGER :: nxl, nxg, ixlg, nxlm1, nxlm2, nxlm3, nxlm4 &
           , nyl, nyg, iylg, nylm1, nylm2, nylm3, nylm4 &
           , nyt,      iytg &
           , nzl &
           , nxlmax, nytmax

  INTEGER, DIMENSION (nprocx) :: nxlp, nylp, ixlgp, iylgp, iprow, ipcol
  INTEGER, DIMENSION (nprocx) :: nytp, iytgp

END MODULE dom_distr_mod

MODULE parallel_mod
  IMPLICIT NONE
  SAVE
  INTEGER :: mpime, nproc, root, group, nprow
  INTEGER :: npcol, myrow, mycol
  LOGICAL :: parallel_build
END MODULE

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
subroutine prtmap(lundia    ,error     ,prsmap    ,lunprt    ,nuprpg    , &
                & nuprln    ,header    ,ipmapc    ,julday    ,dtsec     , &
                & grdang    ,nmax      ,mmax      ,kmax      ,nmaxus    , &
                & lstsci    ,ltur      ,lmaxd     ,lsal      ,ltem      , &
                & namcon    ,kfu       ,kfv       ,kcs       ,s1        , &
                & u1        ,v1        ,wphy      ,alfas     ,r1        , &
                & rtur1     ,vicww     ,dicww     ,rich      ,rho       , &
                & rbuff0    ,rbuff1    ,velt      ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: prtmap.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/prtmap.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Print the time varying map data
!              Selection is done using PRSMAP. For elements like
!              WPHY where KMAX must be > 1 this coupling between
!              KMAX and PRSMAP is done in subroutine RDPRFL
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
!
! Local parameters
!
    integer, parameter :: ncolp = 131
    integer, parameter :: nump = (ncolp - 5)/9
!
! Global variables
!
    integer                                                                    , intent(in)  :: ipmapc !!  Current time counter for printing map
                                                                                !!  data (called with IPMAP (NPMAPC))
    integer                                                                                :: julday !  Description and declaration in inttim.igs
    integer                                                                    , intent(in)  :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                    , intent(in)  :: lmaxd  !  Description and declaration in dimens.igs
    integer                                                                    , intent(in)  :: lsal   !  Description and declaration in dimens.igs
    integer                                                                    , intent(in)  :: lstsci !  Description and declaration in esm_alloc_int.f90
    integer                                                                    , intent(in)  :: ltem   !  Description and declaration in dimens.igs
    integer                                                                    , intent(in)  :: ltur   !  Description and declaration in esm_alloc_int.f90
    integer                                                                                :: lundia !  Description and declaration in inout.igs
    integer                                                                                :: lunprt !  Description and declaration in inout.igs
    integer                                                                    , intent(in)  :: mmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                                :: nmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                    , intent(in)  :: nmaxus !  Description and declaration in esm_alloc_int.f90
    integer                                                                                :: nuprln !  Description and declaration in postpr.igs
    integer                                                                                :: nuprpg !  Description and declaration in postpr.igs
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)                             :: kcs    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               , intent(in)  :: kfu    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               , intent(in)  :: kfv    !  Description and declaration in esm_alloc_int.f90
    logical                                                                                :: error  !!  Flag=TRUE if an error is encountered
    real(fp)                                                                   , intent(in)  :: dtsec  !!  Integration time step [in seconds]
    real(fp)                                                                                :: grdang !  Description and declaration in tricom.igs
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)                            :: alfas  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)              , intent(in)  :: s1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)      , intent(in)  :: dicww  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)      , intent(in)  :: rich   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)      , intent(in)  :: vicww  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax, ltur), intent(in)  :: rtur1  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(in)  :: rho    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(in)  :: u1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(in)  :: v1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(in)  :: wphy   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax, lstsci), intent(in)  :: r1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmaxus, mmax, 0:kmax, lmaxd)                                         :: rbuff0 !!  Help arrays for writing NEFIS files
    real(fp), dimension(nmaxus, mmax, 1:kmax, lmaxd)                                         :: rbuff1 !!  Help arrays for writing NEFIS files
    character(131), dimension(10)                                                            :: header !  Description and declaration in postpr.igs
    character(19)                                                              , intent(in)  :: prsmap !  Description and declaration in tricom.igs
    character(20), dimension(lstsci + ltur)                                    , intent(in)  :: namcon !  Description and declaration in esm_alloc_char.f90
    character(*)                                                               , intent(in)  :: velt   !  Velocity type 'eulerian' or 'GLM'
!
! Local variables
!
    integer       :: idate   ! Absolute date related to ITDATE and TIMSEC
    integer       :: iday     ! Number of days in idate
    integer       :: ihour    ! Number of hours in itime
    integer       :: imin     ! Number of minutes in itime
    integer       :: imon     ! Number of months in idate
    integer       :: isec     ! Number of seconds in itime
    integer       :: itime    ! Absolute time related to ITDATE and TIMSEC
    integer       :: iyear    ! Number of years in idate
    integer       :: k        ! Help var.
    integer       :: kfs      ! Mask for the zeta points (time dependent) =0 inactive point =1 active   point
    integer       :: l        ! Help var.
    integer       :: m        ! Help var.
    integer       :: maxlin   ! Maximum number of lines in one page
    integer       :: md
    integer       :: n        ! Help var.
    integer       :: nd
    integer       :: npr1
    integer       :: npr2
    real(fp)      :: timmin   ! Time since ITDATE in minutes
    real(fp)      :: timsec   ! Time in seconds
    character(20) :: adtsim   ! Absolute date and time related to ITDATE and TIMSEC
    character(8)  :: eenh     ! Unit for concentrations depending on type of constituent
    !
    data maxlin/67/
!
!! executable statements -------------------------------------------------------
!
    ! Define header data 'IPMAPC'
    !
    timsec = ipmapc*dtsec
    timmin = timsec/60.
    !
    call timdat(julday    ,timsec    ,idate    ,itime     )
    iyear = idate/10000
    imon  = (idate - iyear*10000)/100
    iday  = idate - iyear*10000 - imon*100
    ihour = itime/10000
    imin  = (itime - ihour*10000)/100
    isec  = itime - ihour*10000 - imin*100
    write (adtsim, '(i4.4,2(a1,i2.2),a1,3(i2.2,a1))') &
        & iyear, '-', imon, '-', iday, ' ', ihour, ':', imin, ':', isec, ' '
    !
    ! Define header
    !
    write (header(8),'(a,i6,a,f8.2,a,a10,a,a10,a)') &
        & '                                     at NST=',ipmapc,' ( =', timmin, &
        & ' Elapsed min. ,Date : ', adtsim(:10),' time :',adtsim(11:),')'
    !
    ! element 'S1' only if PRSMAP( 1: 1) = 'Y'
    !
    if (prsmap(1:1)=='Y') then
       !
       ! Define water-levels for active zeta points in CM
       !
       do m = 1, mmax
          md = max(1, m - 1)
          do n = 1, nmaxus
             nd = max(1, n - 1)
             kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
             rbuff1(n, m, 1, 1) = s1(n, m)*100.*kfs
          enddo
       enddo
       !
       ! Define text header
       !
       header(7) = ' Computed Water-levels (cm)         '
       !
       ! In sets of NUMP N indices
       !
       do npr1 = 1, nmaxus, nump
          npr2 = min(npr1 + nump - 1, nmaxus)
          !
          ! Test for new page and write HEADER to print file
          !
          write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
          if (nuprln + 6>=maxlin) then
             call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
          else
             write (lunprt, '(a)') header(6)
             write (lunprt, '(a)') header(7)
             write (lunprt, '(a)') header(8)
             write (lunprt, '(a)') header(9)
             write (lunprt, '(a)') header(10)
             nuprln = nuprln + 5
          endif
          !
          ! Print for all M indices
          !
          do m = 1, mmax
             nuprln = nuprln + 1
             !
             ! Test for maximum lines on page
             !
             if (nuprln>maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             endif
             write (lunprt, '(1x,i3,1x,14f9.2)') m, (rbuff1(n, m, 1, 1), n = npr1, npr2)
          enddo
       enddo
    endif
    !
    ! element 'U1' only if PRSMAP( 2: 2) = 'Y'
    !
    if (prsmap(2:2)=='Y') then
       !
       ! Define U-KSI velocities in U-points in CM/SEC
       !
       do k = 1, kmax
          do m = 1, mmax
             do n = 1, nmaxus
                rbuff1(n, m, k, 1) = u1(n, m, k)*100.*kfu(n, m)
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'(3a,i4)') ' Computed U-KSI ', trim(velt), &
              & ' velocities (cm/sec),     layer number', k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.3)') m, (rbuff1(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    endif
    !
    ! element 'V1' only if PRSMAP( 3: 3) = 'Y'
    !
    if (prsmap(3:3)=='Y') then
       !
       ! Define V-ETA velocities in V-points in CM/SEC
       !
       do k = 1, kmax
          do m = 1, mmax
             do n = 1, nmaxus
                rbuff1(n, m, k, 1) = v1(n, m, k)*100.*kfv(n, m)
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'(3a,i4)') ' Computed V-ETA ', trim(velt), &
              & ' velocities (cm/sec),     layer number', k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.3)') m, (rbuff1(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    endif
    !
    ! element 'WPHY' only if KMAX > 1 (:=  PRSMAP( 6: 6) = 'Y')
    !
    if (prsmap(6:6)=='Y') then
       !
       ! Define WPHY velocities in Zeta points in CM/SEC
       !
       do k = 1, kmax
          do m = 1, mmax
             md = max(1, m - 1)
             do n = 1, nmaxus
                nd = max(1, n - 1)
                kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
                rbuff1(n, m, k, 1) = wphy(n, m, k)*100.*kfs
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'('' Computed W-(phys.) velocities (cm/sec), layer number'',i4)') k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.5)') m, (rbuff1(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    endif
    !
    ! element 'R1', only if LSTSCI > 0 (:= PRSMAP( 7:14) <> 'NNNNNNNN')
    !
    if (index(prsmap(7:14), 'Y')/=0) then
       !
       ! Define Concentrations in Zeta-points in variable units
       !
       do l = 1, lstsci
          do k = 1, kmax
             do m = 1, mmax
                md = max(1, m - 1)
                do n = 1, nmaxus
                   nd = max(1, n - 1)
                   kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)&
                       & ))
                   rbuff1(n, m, k, l) = r1(n, m, k, l)*kfs
                enddo
             enddo
          enddo
       enddo
       !
       ! For all constituents, define units
       !
       do l = 1, lstsci
          eenh = '(kg/m3) '
          if (l==lsal) eenh = '(ppt)  '
          if (l==ltem) eenh = '(deg)  '
          !
          ! For all layers
          !
          do k = 1, kmax
             !
             ! Define text header
             !
             write (header(7),'('' Computed '',a,'' '',a,'', layer number'',i4)') &
                 & namcon(l), eenh, k
             !
             ! In sets of NUMP N indices
             !
             do npr1 = 1, nmaxus, nump
                npr2 = min(npr1 + nump - 1, nmaxus)
                !
                ! Test for new page and write HEADER to print file
                !
                write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
                if (nuprln + 6>=maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                else
                   write (lunprt, '(a)') header(6)
                   write (lunprt, '(a)') header(7)
                   write (lunprt, '(a)') header(8)
                   write (lunprt, '(a)') header(9)
                   write (lunprt, '(a)') header(10)
                   nuprln = nuprln + 5
                endif
                !
                ! Print for all M indices
                !
                do m = 1, mmax
                   nuprln = nuprln + 1
                   !
                   ! Test for maximum lines on page
                   !
                   if (nuprln>maxlin) then
                      call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                   endif
                   write (lunprt, '(1x,i3,1x,14f9.3)') m, (rbuff1(n, m, k, l), n = npr1, npr2)
                enddo
             enddo
          enddo
       enddo
    endif
    !
    ! element 'RTUR1', only if LTUR > 0 (:= PRSMAP(15:16) <> 'NN')
    !
    if (index(prsmap(15:16), 'Y')/=0) then
       !
       ! Define Turbulence quantities in Zeta-points
       !
       do l = 1, ltur
          do k = 0, kmax
             do m = 1, mmax
                md = max(1, m - 1)
                do n = 1, nmaxus
                   nd = max(1, n - 1)
                   kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)&
                       & ))
                   rbuff0(n, m, k, l) = rtur1(n, m, k, l)*10000.*kfs
                enddo
             enddo
          enddo
       enddo
       !
       ! For all turbulence quantities
       !
       do l = 1, ltur
          if (l==1) eenh = '(cm2/s3)'
          if (l==2) eenh = '(cm2/s2)'
          !
          ! For all layers
          !
          do k = 0, kmax
             !
             ! Define text header
             !
             write (header(7),'('' Computed '',a,'' '',a,'', layer number'',i4)') &
                 & namcon(lstsci + l), eenh, k
             !
             ! In sets of NUMP N indices
             !
             do npr1 = 1, nmaxus, nump
                npr2 = min(npr1 + nump - 1, nmaxus)
                !
                ! Test for new page and write HEADER to print file
                !
                write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
                if (nuprln + 6>=maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                else
                   write (lunprt, '(a)') header(6)
                   write (lunprt, '(a)') header(7)
                   write (lunprt, '(a)') header(8)
                   write (lunprt, '(a)') header(9)
                   write (lunprt, '(a)') header(10)
                   nuprln = nuprln + 5
                endif
                !
                ! Print for all M indices
                !
                do m = 1, mmax
                   nuprln = nuprln + 1
                   !
                   ! Test for maximum lines on page
                   !
                   if (nuprln>maxlin) then
                      call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                   endif
                   write (lunprt, '(1x,i3,1x,14f9.2)') m, (rbuff0(n, m, k, l), n = npr1, npr2)
                enddo
             enddo
          enddo
       enddo
    endif
    !
    ! element 'VICWW' if KMAX > 1 (:= PRSMAP(17:17) = 'Y')
    !
    if (prsmap(17:17)=='Y') then
       !
       ! Define Vertical eddy Viscosity in Zeta points in CM**2/SEC
       !
       do k = 0, kmax
          do m = 1, mmax
             md = max(1, m - 1)
             do n = 1, nmaxus
                nd = max(1, n - 1)
                kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
                rbuff0(n, m, k, 1) = vicww(n, m, k)*10000.*kfs
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 0, kmax
          !
          ! Define text header
          !
          write (header(7),'('' Computed Vert. Eddy Viscosity  (cm2/s), layer number'',i4)') k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.5)') m, (rbuff0(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    !
    endif
    !
    ! element 'DICWW' if KMAX > 1 (:= PRSMAP(18:18) = 'Y')
    !
    if (prsmap(18:18)=='Y') then
       !
       ! Define Vertical Eddy Diffusivity in Zeta points in CM**2/SEC
       !
       do k = 0, kmax
          do m = 1, mmax
             md = max(1, m - 1)
             do n = 1, nmaxus
                nd = max(1, n - 1)
                kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
                rbuff0(n, m, k, 1) = dicww(n, m, k)*10000.*kfs
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'('' Computed Vert.Eddy Diffusivity (cm2/s), layer number'',i4)') k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.5)') m, (rbuff0(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    !
    endif
    !
    ! element 'RICH' if KMAX > 1 (:= PRSMAP(17:18) <> 'NN')
    !
    if (index(prsmap(17:18), 'Y')>0) then
       !
       ! Define Richardson Numbers in zeta points
       !
       do k = 0, kmax
          do m = 1, mmax
             md = max(1, m - 1)
             do n = 1, nmaxus
                nd = max(1, n - 1)
                kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
                rbuff0(n, m, k, 1) = rich(n, m, k)*kfs
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'(a,i4)') ' Computed Richardson numbers(-), layer number', k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.5)') m, (rbuff0(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    !
    endif
    !
    ! element 'RHO' if LSAL > 0 or LTEM > 0 (:= PRSMAP(19:19) = 'Y')
    !
    if (prsmap(19:19)=='Y') then
       !
       ! Define Densities in Zeta points in KG/M**3
       !
       do k = 1, kmax
          do m = 1, mmax
             md = max(1, m - 1)
             do n = 1, nmaxus
                nd = max(1, n - 1)
                kfs = min(1, (kfu(n, m) + kfu(n, md) + kfv(n, m) + kfv(nd, m)))
                rbuff1(n, m, k, 1) = rho(n, m, k)*kfs
             enddo
          enddo
       enddo
       !
       ! For all layers
       !
       do k = 1, kmax
          !
          ! Define text header
          !
          write (header(7),'('' Computed Densities             (kg/m3), layer number'',i4)') k
          !
          ! In sets of NUMP N indices
          !
          do npr1 = 1, nmaxus, nump
             npr2 = min(npr1 + nump - 1, nmaxus)
             !
             ! Test for new page and write HEADER to print file
             !
             write (header(9), '(''   N='',14i9)') (n, n = npr1, npr2)
             if (nuprln + 6>=maxlin) then
                call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
             else
                write (lunprt, '(a)') header(6)
                write (lunprt, '(a)') header(7)
                write (lunprt, '(a)') header(8)
                write (lunprt, '(a)') header(9)
                write (lunprt, '(a)') header(10)
                nuprln = nuprln + 5
             endif
             !
             ! Print for all M indices
             !
             do m = 1, mmax
                nuprln = nuprln + 1
                !
                ! Test for maximum lines on page
                !
                if (nuprln>maxlin) then
                   call prhead(lunprt    ,header    ,nuprpg    ,nuprln    )
                endif
                write (lunprt, '(1x,i3,1x,14f9.2)') m, (rbuff1(n, m, k, 1), n = npr1, npr2)
             enddo
          enddo
       enddo
    !
    endif
end subroutine prtmap

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd