solver.F90 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133
  1. MODULE solver
  2. !!======================================================================
  3. !! *** MODULE solver ***
  4. !! Ocean solver : initialization of ocean solver
  5. !!=====================================================================
  6. !! History : OPA ! 1990-10 (G. Madec) Original code
  7. !! ! 1993-02 (O. Marti)
  8. !! ! 1997-02 (G. Madec) local depth inverse computation
  9. !! ! 1998-10 (G. Roullet, G. Madec) free surface
  10. !! NEMO 1.0 ! 2003-07 (G. Madec) free form, F90
  11. !! 3.2 ! 2009-07 (R. Benshila) suppression of rigid-lid & FETI solver
  12. !!----------------------------------------------------------------------
  13. #if defined key_dynspg_flt || defined key_esopa
  14. !!----------------------------------------------------------------------
  15. !! 'key_dynspg_flt' filtered free surface
  16. !!----------------------------------------------------------------------
  17. !! solver_init: solver initialization
  18. !!----------------------------------------------------------------------
  19. USE oce ! ocean dynamics and tracers variables
  20. USE dom_oce ! ocean space and time domain variables
  21. USE zdf_oce ! ocean vertical physics variables
  22. USE sol_oce ! solver variables
  23. USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
  24. USE solmat ! matrix of the solver
  25. USE in_out_manager ! I/O manager
  26. USE lbclnk ! ocean lateral boundary conditions (or mpp link)
  27. USE lib_mpp ! MPP library
  28. USE timing ! timing
  29. IMPLICIT NONE
  30. !!----------------------------------------------------------------------
  31. !! NEMO/OPA 3.3 , NEMO Consortium (2010)
  32. !! $Id: solver.F90 5533 2015-07-02 13:50:50Z jchanut $
  33. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  34. !!----------------------------------------------------------------------
  35. CONTAINS
  36. SUBROUTINE solver_init( kt )
  37. !!----------------------------------------------------------------------
  38. !! *** ROUTINE solver_init ***
  39. !!
  40. !! ** Purpose : Initialization of the elliptic solver
  41. !!
  42. !! ** Method : a solver is required when using the filtered free
  43. !! surface.
  44. !!
  45. !! ** Action : - c_solver_pt : nature of the gridpoint at which the solver is applied
  46. !!
  47. !! References : Jensen, 1986: Adv. Phys. Oceanogr. Num. Mod.,Ed. O Brien,87-110.
  48. !!----------------------------------------------------------------------
  49. INTEGER, INTENT(in) :: kt
  50. INTEGER :: ios ! Local integer output status for namelist read
  51. !!
  52. NAMELIST/namsol/ nn_solv, nn_sol_arp, nn_nmin, nn_nmax, nn_nmod, rn_eps, rn_resmax, rn_sor
  53. !!----------------------------------------------------------------------
  54. !
  55. IF( nn_timing == 1 ) CALL timing_start('solver_init')
  56. !
  57. IF(lwp) THEN !* open elliptic solver statistics file (only on the printing processors)
  58. CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
  59. ENDIF
  60. REWIND( numnam_ref ) ! Namelist namsol in reference namelist : Elliptic solver / free surface
  61. READ ( numnam_ref, namsol, IOSTAT = ios, ERR = 901)
  62. 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsol in reference namelist', lwp )
  63. REWIND( numnam_cfg ) ! Namelist namsol in configuration namelist : Elliptic solver / free surface
  64. READ ( numnam_cfg, namsol, IOSTAT = ios, ERR = 902 )
  65. 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsol in configuration namelist', lwp )
  66. IF(lwm) WRITE ( numond, namsol )
  67. IF(lwp) THEN !* Control print
  68. WRITE(numout,*)
  69. WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient'
  70. WRITE(numout,*) '~~~~~~~~~~~'
  71. WRITE(numout,*) ' Namelist namsol : set solver parameters'
  72. WRITE(numout,*) ' type of elliptic solver nn_solv = ', nn_solv
  73. WRITE(numout,*) ' absolute/relative (0/1) precision nn_sol_arp = ', nn_sol_arp
  74. WRITE(numout,*) ' minimum iterations for solver nn_nmin = ', nn_nmin
  75. WRITE(numout,*) ' maximum iterations for solver nn_nmax = ', nn_nmax
  76. WRITE(numout,*) ' frequency for test nn_nmod = ', nn_nmod
  77. WRITE(numout,*) ' absolute precision of solver rn_eps = ', rn_eps
  78. WRITE(numout,*) ' absolute precision for SOR solver rn_resmax = ', rn_resmax
  79. WRITE(numout,*) ' optimal coefficient of sor rn_sor = ', rn_sor
  80. WRITE(numout,*)
  81. ENDIF
  82. eps = rn_eps
  83. ! ! allocate solver arrays
  84. IF( .NOT. lk_agrif .OR. .NOT. ln_rstart) THEN
  85. IF( sol_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'solver_init : unable to allocate sol_oce arrays' )
  86. gcx (:,:) = 0.e0
  87. gcxb(:,:) = 0.e0
  88. ENDIF
  89. SELECT CASE( nn_solv ) !* parameter check
  90. !
  91. CASE ( 1 ) ! preconditioned conjugate gradient solver
  92. IF(lwp) WRITE(numout,*) ' a preconditioned conjugate gradient solver is used'
  93. IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) CALL ctl_stop( ' jpr2di and jpr2dj should be equal to zero' )
  94. !
  95. CASE ( 2 ) ! successive-over-relaxation solver
  96. IF(lwp) WRITE(numout,*) ' a successive-over-relaxation solver with extra outer halo is used'
  97. IF(lwp) WRITE(numout,*) ' with jpr2di =', jpr2di, ' and jpr2dj =', jpr2dj
  98. IF( .NOT. lk_mpp .AND. jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
  99. CALL ctl_stop( 'jpr2di and jpr2dj are not equal to zero', &
  100. & 'In this case the algorithm should be used only with the key_mpp_... option' )
  101. ELSE
  102. IF( ( ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) .OR. ( jpni /= 1 ) ) &
  103. & .AND. ( jpr2di /= jpr2dj ) ) CALL ctl_stop( 'jpr2di should be equal to jpr2dj' )
  104. ENDIF
  105. !
  106. CASE DEFAULT ! error in parameter
  107. WRITE(ctmp1,*) ' bad flag value for nn_solv = ', nn_solv
  108. CALL ctl_stop( ctmp1 )
  109. END SELECT
  110. !
  111. ! !* Grid-point at which the solver is applied
  112. !!gm c_solver_pt should be removed: nomore bsf, only T-point is used
  113. c_solver_pt = 'T' ! always T-point (ssh solver only, not anymore bsf)
  114. CALL sol_mat( kt ) !* Construction of the elliptic system matrix
  115. !
  116. IF( nn_timing == 1 ) CALL timing_stop('solver_init')
  117. !
  118. END SUBROUTINE solver_init
  119. #endif
  120. !!======================================================================
  121. END MODULE solver