diaharm.F90 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565
  1. MODULE diaharm
  2. !!======================================================================
  3. !! *** MODULE diaharm ***
  4. !! Harmonic analysis of tidal constituents
  5. !!======================================================================
  6. !! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code
  7. !!----------------------------------------------------------------------
  8. #if defined key_diaharm && defined key_tide
  9. !!----------------------------------------------------------------------
  10. !! 'key_diaharm'
  11. !! 'key_tide'
  12. !!----------------------------------------------------------------------
  13. USE oce ! ocean dynamics and tracers variables
  14. USE dom_oce ! ocean space and time domain
  15. USE phycst
  16. USE dynspg_oce
  17. USE dynspg_ts
  18. USE daymod
  19. USE tide_mod
  20. !
  21. USE in_out_manager ! I/O units
  22. USE iom ! I/0 library
  23. USE ioipsl ! NetCDF IPSL library
  24. USE lbclnk ! ocean lateral boundary conditions (or mpp link)
  25. USE diadimg ! To write dimg
  26. USE timing ! preformance summary
  27. USE wrk_nemo ! working arrays
  28. IMPLICIT NONE
  29. PRIVATE
  30. LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
  31. INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
  32. INTEGER, PARAMETER :: jpdimsparse = jpincomax*300*24
  33. ! !!** namelist variables **
  34. INTEGER :: nit000_han ! First time step used for harmonic analysis
  35. INTEGER :: nitend_han ! Last time step used for harmonic analysis
  36. INTEGER :: nstep_han ! Time step frequency for harmonic analysis
  37. INTEGER :: nb_ana ! Number of harmonics to analyse
  38. INTEGER , ALLOCATABLE, DIMENSION(:) :: name
  39. REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
  40. REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
  41. REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
  42. INTEGER :: ninco, nsparse
  43. INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
  44. INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
  45. REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
  46. REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
  47. REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
  48. REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
  49. CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
  50. PUBLIC dia_harm ! routine called by step.F90
  51. !!----------------------------------------------------------------------
  52. !! NEMO/OPA 3.5 , NEMO Consortium (2013)
  53. !! $Id: diaharm.F90 2544 2015-08-24 09:00:45Z ufla $
  54. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  55. !!----------------------------------------------------------------------
  56. CONTAINS
  57. SUBROUTINE dia_harm_init
  58. !!----------------------------------------------------------------------
  59. !! *** ROUTINE dia_harm_init ***
  60. !!
  61. !! ** Purpose : Initialization of tidal harmonic analysis
  62. !!
  63. !! ** Method : Initialize frequency array and nodal factor for nit000_han
  64. !!
  65. !!--------------------------------------------------------------------
  66. INTEGER :: jh, nhan, jk, ji
  67. INTEGER :: ios ! Local integer output status for namelist read
  68. NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
  69. !!----------------------------------------------------------------------
  70. IF(lwp) THEN
  71. WRITE(numout,*)
  72. WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
  73. WRITE(numout,*) '~~~~~~~ '
  74. ENDIF
  75. !
  76. CALL tide_init_Wave
  77. !
  78. REWIND( numnam_ref ) ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
  79. READ ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
  80. 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
  81. REWIND( numnam_cfg ) ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
  82. READ ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
  83. 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
  84. IF(lwm) WRITE ( numond, nam_diaharm )
  85. !
  86. IF(lwp) THEN
  87. WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han
  88. WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han
  89. WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han
  90. ENDIF
  91. ! Basic checks on harmonic analysis time window:
  92. ! ----------------------------------------------
  93. IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', &
  94. & ' restart capability not implemented' )
  95. IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', &
  96. & 'restart capability not implemented' )
  97. IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) &
  98. & CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
  99. nb_ana = 0
  100. DO jk=1,jpmax_harmo
  101. DO ji=1,jpmax_harmo
  102. IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
  103. nb_ana=nb_ana+1
  104. ENDIF
  105. END DO
  106. END DO
  107. !
  108. IF(lwp) THEN
  109. WRITE(numout,*) ' Namelist nam_diaharm'
  110. WRITE(numout,*) ' nb_ana = ', nb_ana
  111. CALL flush(numout)
  112. ENDIF
  113. !
  114. IF (nb_ana > jpmax_harmo) THEN
  115. IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
  116. IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
  117. nstop = nstop + 1
  118. ENDIF
  119. ALLOCATE(name (nb_ana))
  120. DO jk=1,nb_ana
  121. DO ji=1,jpmax_harmo
  122. IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
  123. name(jk) = ji
  124. EXIT
  125. END IF
  126. END DO
  127. END DO
  128. ! Initialize frequency array:
  129. ! ---------------------------
  130. ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
  131. CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
  132. IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
  133. DO jh = 1, nb_ana
  134. IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh)
  135. END DO
  136. ! Initialize temporary arrays:
  137. ! ----------------------------
  138. ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
  139. ana_temp(:,:,:,:) = 0._wp
  140. END SUBROUTINE dia_harm_init
  141. SUBROUTINE dia_harm ( kt )
  142. !!----------------------------------------------------------------------
  143. !! *** ROUTINE dia_harm ***
  144. !!
  145. !! ** Purpose : Tidal harmonic analysis main routine
  146. !!
  147. !! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han]
  148. !!
  149. !!--------------------------------------------------------------------
  150. INTEGER, INTENT( IN ) :: kt
  151. !
  152. INTEGER :: ji, jj, jh, jc, nhc
  153. REAL(wp) :: ztime, ztemp
  154. !!--------------------------------------------------------------------
  155. IF( nn_timing == 1 ) CALL timing_start('dia_harm')
  156. IF( kt == nit000 ) CALL dia_harm_init
  157. IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
  158. ztime = (kt-nit000+1) * rdt
  159. nhc = 0
  160. DO jh = 1, nb_ana
  161. DO jc = 1, 2
  162. nhc = nhc+1
  163. ztemp =( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
  164. & +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
  165. DO jj = 1,jpj
  166. DO ji = 1,jpi
  167. ! Elevation
  168. ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)
  169. ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj)
  170. ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj)
  171. END DO
  172. END DO
  173. !
  174. END DO
  175. END DO
  176. !
  177. END IF
  178. IF ( kt == nitend_han ) CALL dia_harm_end
  179. IF( nn_timing == 1 ) CALL timing_stop('dia_harm')
  180. END SUBROUTINE dia_harm
  181. SUBROUTINE dia_harm_end
  182. !!----------------------------------------------------------------------
  183. !! *** ROUTINE diaharm_end ***
  184. !!
  185. !! ** Purpose : Compute the Real and Imaginary part of tidal constituents
  186. !!
  187. !! ** Action : Decompose the signal on the harmonic constituents
  188. !!
  189. !!--------------------------------------------------------------------
  190. INTEGER :: ji, jj, jh, jc, jn, nhan, jl
  191. INTEGER :: ksp, kun, keq
  192. REAL(wp) :: ztime, ztime_ini, ztime_end
  193. REAL(wp) :: X1,X2
  194. REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
  195. !!--------------------------------------------------------------------
  196. CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
  197. IF(lwp) WRITE(numout,*)
  198. IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
  199. IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
  200. ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
  201. ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
  202. nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
  203. ninco = 2*nb_ana
  204. ksp = 0
  205. keq = 0
  206. DO jn = 1, nhan
  207. ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
  208. keq = keq + 1
  209. kun = 0
  210. DO jh = 1, nb_ana
  211. DO jc = 1, 2
  212. kun = kun + 1
  213. ksp = ksp + 1
  214. nisparse(ksp) = keq
  215. njsparse(ksp) = kun
  216. valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
  217. & + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
  218. END DO
  219. END DO
  220. END DO
  221. nsparse = ksp
  222. ! Elevation:
  223. DO jj = 1, jpj
  224. DO ji = 1, jpi
  225. ! Fill input array
  226. kun = 0
  227. DO jh = 1, nb_ana
  228. DO jc = 1, 2
  229. kun = kun + 1
  230. ztmp4(kun)=ana_temp(ji,jj,kun,1)
  231. END DO
  232. END DO
  233. CALL SUR_DETERMINE(jj)
  234. ! Fill output array
  235. DO jh = 1, nb_ana
  236. ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
  237. ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
  238. END DO
  239. END DO
  240. END DO
  241. ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
  242. & out_u (jpi,jpj,2*nb_ana), &
  243. & out_v (jpi,jpj,2*nb_ana) )
  244. DO jj = 1, jpj
  245. DO ji = 1, jpi
  246. DO jh = 1, nb_ana
  247. X1 = ana_amp(ji,jj,jh,1)
  248. X2 =-ana_amp(ji,jj,jh,2)
  249. out_eta(ji,jj,jh ) = X1 * tmask_i(ji,jj)
  250. out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
  251. END DO
  252. END DO
  253. END DO
  254. ! ubar:
  255. DO jj = 1, jpj
  256. DO ji = 1, jpi
  257. ! Fill input array
  258. kun=0
  259. DO jh = 1,nb_ana
  260. DO jc = 1,2
  261. kun = kun + 1
  262. ztmp4(kun)=ana_temp(ji,jj,kun,2)
  263. END DO
  264. END DO
  265. CALL SUR_DETERMINE(jj+1)
  266. ! Fill output array
  267. DO jh = 1, nb_ana
  268. ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
  269. ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
  270. END DO
  271. END DO
  272. END DO
  273. DO jj = 1, jpj
  274. DO ji = 1, jpi
  275. DO jh = 1, nb_ana
  276. X1= ana_amp(ji,jj,jh,1)
  277. X2=-ana_amp(ji,jj,jh,2)
  278. out_u(ji,jj, jh) = X1 * umask_i(ji,jj)
  279. out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
  280. ENDDO
  281. ENDDO
  282. ENDDO
  283. ! vbar:
  284. DO jj = 1, jpj
  285. DO ji = 1, jpi
  286. ! Fill input array
  287. kun=0
  288. DO jh = 1,nb_ana
  289. DO jc = 1,2
  290. kun = kun + 1
  291. ztmp4(kun)=ana_temp(ji,jj,kun,3)
  292. END DO
  293. END DO
  294. CALL SUR_DETERMINE(jj+1)
  295. ! Fill output array
  296. DO jh = 1, nb_ana
  297. ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
  298. ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
  299. END DO
  300. END DO
  301. END DO
  302. DO jj = 1, jpj
  303. DO ji = 1, jpi
  304. DO jh = 1, nb_ana
  305. X1=ana_amp(ji,jj,jh,1)
  306. X2=-ana_amp(ji,jj,jh,2)
  307. out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)
  308. out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
  309. END DO
  310. END DO
  311. END DO
  312. CALL dia_wri_harm ! Write results in files
  313. CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
  314. !
  315. END SUBROUTINE dia_harm_end
  316. SUBROUTINE dia_wri_harm
  317. !!--------------------------------------------------------------------
  318. !! *** ROUTINE dia_wri_harm ***
  319. !!
  320. !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
  321. !!--------------------------------------------------------------------
  322. CHARACTER(LEN=lc) :: cltext
  323. CHARACTER(LEN=lc) :: &
  324. cdfile_name_T , & ! name of the file created (T-points)
  325. cdfile_name_U , & ! name of the file created (U-points)
  326. cdfile_name_V ! name of the file created (V-points)
  327. INTEGER :: jh
  328. !!----------------------------------------------------------------------
  329. #if defined key_dimgout
  330. cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
  331. cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
  332. cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
  333. #endif
  334. IF(lwp) WRITE(numout,*) ' '
  335. IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
  336. #if defined key_dimgout
  337. IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)
  338. IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)
  339. IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)
  340. #endif
  341. IF(lwp) WRITE(numout,*) ' '
  342. ! A) Elevation
  343. !/////////////
  344. !
  345. #if defined key_dimgout
  346. cltext='Elevation amplitude and phase'
  347. CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
  348. #else
  349. DO jh = 1, nb_ana
  350. CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
  351. CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
  352. END DO
  353. #endif
  354. ! B) ubar
  355. !/////////
  356. !
  357. #if defined key_dimgout
  358. cltext='ubar amplitude and phase'
  359. CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
  360. #else
  361. DO jh = 1, nb_ana
  362. CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
  363. CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
  364. END DO
  365. #endif
  366. ! C) vbar
  367. !/////////
  368. !
  369. #if defined key_dimgout
  370. cltext='vbar amplitude and phase'
  371. CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
  372. #else
  373. DO jh = 1, nb_ana
  374. CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) )
  375. CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
  376. END DO
  377. #endif
  378. !
  379. END SUBROUTINE dia_wri_harm
  380. SUBROUTINE SUR_DETERMINE(init)
  381. !!---------------------------------------------------------------------------------
  382. !! *** ROUTINE SUR_DETERMINE ***
  383. !!
  384. !!
  385. !!
  386. !!---------------------------------------------------------------------------------
  387. INTEGER, INTENT(in) :: init
  388. !
  389. INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
  390. REAL(wp) :: zval1, zval2, zx1
  391. REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
  392. INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
  393. !---------------------------------------------------------------------------------
  394. CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
  395. CALL wrk_alloc( jpincomax , ipos2 , ipivot )
  396. IF( init == 1 ) THEN
  397. IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
  398. IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
  399. !
  400. ztmp3(:,:) = 0._wp
  401. !
  402. DO jk1_sd = 1, nsparse
  403. DO jk2_sd = 1, nsparse
  404. nisparse(jk2_sd) = nisparse(jk2_sd)
  405. njsparse(jk2_sd) = njsparse(jk2_sd)
  406. IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
  407. ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
  408. & + valuesparse(jk1_sd)*valuesparse(jk2_sd)
  409. ENDIF
  410. END DO
  411. END DO
  412. !
  413. DO jj_sd = 1 ,ninco
  414. ipos1(jj_sd) = jj_sd
  415. ipos2(jj_sd) = jj_sd
  416. END DO
  417. !
  418. DO ji_sd = 1 , ninco
  419. !
  420. !find greatest non-zero pivot:
  421. zval1 = ABS(ztmp3(ji_sd,ji_sd))
  422. !
  423. ipivot(ji_sd) = ji_sd
  424. DO jj_sd = ji_sd, ninco
  425. zval2 = ABS(ztmp3(ji_sd,jj_sd))
  426. IF( zval2.GE.zval1 )THEN
  427. ipivot(ji_sd) = jj_sd
  428. zval1 = zval2
  429. ENDIF
  430. END DO
  431. !
  432. DO ji1_sd = 1, ninco
  433. zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
  434. zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
  435. ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
  436. ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
  437. END DO
  438. !
  439. ipos2(ji_sd) = ipos1(ipivot(ji_sd))
  440. ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
  441. ipos1(ji_sd) = ipos2(ji_sd)
  442. ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
  443. zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
  444. DO jj_sd = 1, ninco
  445. ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
  446. END DO
  447. !
  448. DO ji2_sd = ji_sd+1, ninco
  449. zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
  450. DO jj_sd=1,ninco
  451. ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
  452. END DO
  453. END DO
  454. !
  455. END DO
  456. !
  457. ENDIF ! End init==1
  458. DO ji_sd = 1, ninco
  459. ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
  460. DO ji2_sd = ji_sd+1, ninco
  461. ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
  462. END DO
  463. END DO
  464. !system solving:
  465. ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
  466. ji_sd = ninco
  467. DO ji_sd = ninco-1, 1, -1
  468. zx1 = 0._wp
  469. DO jj_sd = ji_sd+1, ninco
  470. zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
  471. END DO
  472. ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
  473. END DO
  474. DO jj_sd =1, ninco
  475. ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
  476. END DO
  477. CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
  478. CALL wrk_dealloc( jpincomax , ipos2 , ipivot )
  479. !
  480. END SUBROUTINE SUR_DETERMINE
  481. #else
  482. !!----------------------------------------------------------------------
  483. !! Default case : Empty module
  484. !!----------------------------------------------------------------------
  485. LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
  486. CONTAINS
  487. SUBROUTINE dia_harm ( kt ) ! Empty routine
  488. INTEGER, INTENT( IN ) :: kt
  489. WRITE(*,*) 'dia_harm: you should not have seen this print'
  490. END SUBROUTINE dia_harm
  491. #endif
  492. !!======================================================================
  493. END MODULE diaharm