Changeset 1556 for trunk/NEMO/OPA_SRC/SOL/solmat.F90
- Timestamp:
- 2009-07-29T12:51:45+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SOL/solmat.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SOL/solmat.F90
r1528 r1556 9 9 !! ! 96-05 (G. Madec) merge sor and pcg formulations 10 10 !! ! 96-11 (A. Weaver) correction to preconditioning 11 !! ! 98-02 (M. Guyon) FETI method12 11 !! 8.5 ! 02-08 (G. Madec) F90: Free form 13 12 !! ! 02-11 (C. Talandier, A-M. Treguier) Free surface & Open boundaries … … 51 50 !! 52 51 !! ** Purpose : Construction of the matrix of used by the elliptic 53 !! solvers (either sor , pcg or fetimethods).52 !! solvers (either sor or pcg methods). 54 53 !! 55 54 !! ** Method : … … 57 56 !! The matrix is built for the divergence of the transport system 58 57 !! a diagonal preconditioning matrix is also defined. 59 !! Note that for feti solver (nsolv=3) a specific initialization60 !! is required (call to fetstr.F) for memory allocation and inter-61 !! face definition.62 58 !! 63 59 !! ** Action : - gcp : extra-diagonal elements of the matrix … … 74 70 REAL(wp) :: z2dt, zcoef 75 71 !!---------------------------------------------------------------------- 76 77 ! FETI method ( nsolv = 3)78 ! memory allocation and interface definition for the solver79 80 IF( nsolv == 3 ) CALL fetstr81 72 82 73 … … 275 266 ! ------------------------ 276 267 277 IF( nsolv /= 3 ) THEN 268 ! SOR and PCG solvers 269 DO jj = 1, jpj 270 DO ji = 1, jpi 271 IF( bmask(ji,jj) /= 0.e0 ) gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj) 272 END DO 273 END DO 278 274 279 ! SOR and PCG solvers 280 DO jj = 1, jpj 281 DO ji = 1, jpi 282 IF( bmask(ji,jj) /= 0.e0 ) gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj) 283 END DO 284 END DO 285 286 gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:) 287 gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:) 288 gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:) 289 gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:) 290 IF( nsolv == 2 ) gccd(:,:) = sor * gcp(:,:,2) 291 292 IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN 293 CALL lbc_lnk_e( gcp (:,:,1), c_solver_pt, 1. ) ! lateral boundary conditions 294 CALL lbc_lnk_e( gcp (:,:,2), c_solver_pt, 1. ) ! lateral boundary conditions 295 CALL lbc_lnk_e( gcp (:,:,3), c_solver_pt, 1. ) ! lateral boundary conditions 296 CALL lbc_lnk_e( gcp (:,:,4), c_solver_pt, 1. ) ! lateral boundary conditions 297 CALL lbc_lnk_e( gcdprc(:,:) , c_solver_pt, 1. ) ! lateral boundary conditions 298 CALL lbc_lnk_e( gcdmat(:,:) , c_solver_pt, 1. ) ! lateral boundary conditions 299 IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements 300 END IF 301 302 ELSE 303 304 ! FETI method 305 ! if feti solver : gcdprc is a mask for the non-overlapping 306 ! data structuring 307 308 DO jj = 1, jpj 309 DO ji = 1, jpi 310 IF( bmask(ji,jj) /= 0.e0 ) THEN 311 gcdprc(ji,jj) = 1.e0 312 ELSE 313 gcdprc(ji,jj) = 0.e0 314 ENDIF 315 END DO 316 END DO 317 318 ! so "common" line & "common" column have to be !=0 except on global 319 ! domain boundaries 320 ! pbs with nbondi if nperio != 2 ? 321 ! ii = nldi-1 322 ! pb with nldi value if jperio==1 : nbondi modifyed at the end 323 ! of inimpp.F => pb 324 ! pb with periodicity conditions : iiend, ijend 325 326 ijend = nlej 327 iiend = nlei 328 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) iiend = nlci - jpreci 329 ii = jpreci 330 331 ! case number 1 332 333 IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN 334 DO jj = 1, ijend 335 IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1. 336 END DO 337 ENDIF 338 339 ! case number 2 340 341 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 342 DO jj = 1, ijend 343 IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1. 344 END DO 345 ENDIF 346 347 ! ij = nldj-1 348 ! pb with nldi value if jperio==1 : nbondi modifyed at the end 349 ! of inimpp.F => pb, here homogeneisation... 350 351 ij = jprecj 352 IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN 353 DO ji = 1, iiend 354 IF( fmask(ji,ij,1) == 1. ) gcdprc(ji,ij) = 1. 355 END DO 356 ENDIF 357 ENDIF 358 359 275 gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:) 276 gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:) 277 gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:) 278 gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:) 279 IF( nsolv == 2 ) gccd(:,:) = sor * gcp(:,:,2) 280 281 IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN 282 CALL lbc_lnk_e( gcp (:,:,1), c_solver_pt, 1. ) ! lateral boundary conditions 283 CALL lbc_lnk_e( gcp (:,:,2), c_solver_pt, 1. ) ! lateral boundary conditions 284 CALL lbc_lnk_e( gcp (:,:,3), c_solver_pt, 1. ) ! lateral boundary conditions 285 CALL lbc_lnk_e( gcp (:,:,4), c_solver_pt, 1. ) ! lateral boundary conditions 286 CALL lbc_lnk_e( gcdprc(:,:) , c_solver_pt, 1. ) ! lateral boundary conditions 287 CALL lbc_lnk_e( gcdmat(:,:) , c_solver_pt, 1. ) ! lateral boundary conditions 288 IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements 289 END IF 290 360 291 ! 4. Initialization the arrays used in pcg 361 292 ! ---------------------------------------- … … 364 295 gcdes(:,:) = 0.e0 365 296 gccd (:,:) = 0.e0 366 367 ! FETI method 368 IF( nsolv == 3 ) THEN 369 CALL fetmat ! Matrix treatment : Neumann condition, inverse computation 370 CALL fetsch ! data framework for the Schur Dual solver 371 ENDIF 372 297 ! 373 298 END SUBROUTINE sol_mat 374 299 … … 470 395 471 396 END SELECT ! npolj 472 397 ! 473 398 END SUBROUTINE sol_exd 474 475 #if defined key_feti476 477 SUBROUTINE fetstr478 !!---------------------------------------------------------------------479 !! *** ROUTINE fetstr ***480 !!481 !! ** Purpose : Construction of the matrix of the barotropic stream482 !! function system.483 !! Finite Elements Tearing & Interconnecting (FETI) approach484 !! Memory allocation and interface definition for the solver485 !!486 !! ** Method :487 !!488 !! References :489 !! Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :490 !! A domain decomposition solver to compute the barotropic491 !! component of an OGCM in the parallel processing field.492 !! Ocean Modelling, issue 105, december 94.493 !!494 !! History :495 !! ! 98-02 (M. Guyon) Original code496 !! 8.5 ! 02-09 (G. Madec) F90: Free form and module497 !!----------------------------------------------------------------------498 !! * Modules used499 USE lib_feti ! feti librairy500 !! * Local declarations501 INTEGER :: iiend, ijend, iperio ! temporary integers502 !!---------------------------------------------------------------------503 504 505 ! Preconditioning technics of the Dual Schur Operator506 ! <= definition of the Coarse Grid solver507 ! <= dimension of the nullspace of the local operators508 ! <= Neumann boundaries conditions509 510 ! 0. Initializations511 ! ------------------512 513 ndkerep = 1514 515 ! initialization of the superstructures management516 517 malxm = 1518 malim = 1519 520 ! memory space for the pcpg associated with the FETI dual formulation521 ! ndkerep is associated to the list of rigid modes,522 ! ndkerep == 1 because the Dual Operator523 ! is a first order operator due to SPG elliptic Operator is a524 ! second order operator525 526 nim = 50527 nim = nim + ndkerep528 nim = nim + 2*jpi + 2*jpj529 nim = nim + jpi*jpj530 531 nxm = 33532 nxm = nxm + 4*jpnij533 nxm = nxm + 19*(jpi+jpj)534 nxm = nxm + 13*jpi*jpj535 nxm = nxm + jpi*jpi*jpj536 537 ! krylov space memory538 539 iperio = 0540 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6) iperio = 1541 nxm = nxm + 3*(jpnij-jpni)*jpi542 nxm = nxm + 3*(jpnij-jpnj+iperio)*jpj543 nxm = nxm + 2*(jpi+jpj)*(jpnij-jpni)*jpi544 nxm = nxm + 2*(jpi+jpj)*(jpnij-jpnj+iperio)*jpj545 546 ! Resolution with the Schur dual Method ( frontal and local solver by547 ! blocks548 ! Case with a local symetrical matrix549 ! The local matrix is stored in a multi-column form550 ! The total number of nodes for this subdomain is named "noeuds"551 552 noeuds = jpi*jpj553 nifmat = jpi-1554 njfmat = jpj-1555 nelem = nifmat*njfmat556 npe = 4557 nmorse = 5*noeuds558 559 ! 1. mesh building560 ! ----------------561 562 ! definition of specific information for a subdomain563 ! narea : subdomain number = processor number +1564 ! ninterf : neighbour subdomain number565 ! nni : interface point number566 ! ndvois array : neighbour subdomain list567 ! maplistin array : node pointer at each interface568 ! maplistin array : concatened list of interface nodes569 570 ! messag coding is necessary by interface type for avoid collision571 ! if nperio == 1572 573 ! lint array : indoor interface list / type574 ! lext array : outdoor interface list / type575 576 ! domain with jpniXjpnj subdomains577 578 CALL feti_inisub(nifmat,njfmat,nbondi,nbondj,nperio, &579 nbsw,nbnw,nbse,nbne,ninterf,ninterfc,nni,nnic)580 581 CALL feti_creadr(malim,malimax,nim,3*ninterf ,mandvois ,'ndvois' )582 CALL feti_creadr(malim,malimax,nim,3*ninterfc,mandvoisc,'ndvoisc')583 CALL feti_creadr(malim,malimax,nim,ninterfc+1,maplistin,'plistin')584 CALL feti_creadr(malim,malimax,nim,nnic ,malistin ,'listin' )585 586 ! pb with periodicity conditions : iiend, ijend587 588 ijend = nlej589 iiend = nlei590 IF (jperio == 1) iiend = nlci - jpreci591 592 CALL feti_subound(nifmat,njfmat,nldi,iiend,nldj,ijend, &593 narea,nbondi,nbondj,nperio, &594 ninterf,ninterfc, &595 nowe,noea,noso,nono, &596 nbsw,nbnw,nbse,nbne, &597 npsw,npnw,npse,npne, &598 mfet(mandvois),mfet(mandvoisc), &599 mfet(maplistin),nnic,mfet(malistin) )600 601 END SUBROUTINE fetstr602 603 604 SUBROUTINE fetmat605 !!---------------------------------------------------------------------606 !! *** ROUTINE fetmat ***607 !!608 !! ** Purpose : Construction of the matrix of the barotropic stream609 !! function system.610 !! Finite Elements Tearing & Interconnecting (FETI) approach611 !! Matrix treatment : Neumann condition, inverse computation612 !!613 !! ** Method :614 !!615 !! References :616 !! Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :617 !! A domain decomposition solver to compute the barotropic618 !! component of an OGCM in the parallel processing field.619 !! Ocean Modelling, issue 105, december 94.620 !!621 !! History :622 !! ! 98-02 (M. Guyon) Original code623 !! 8.5 ! 02-09 (G. Madec) F90: Free form and module624 !!----------------------------------------------------------------------625 !! * Modules used626 USE lib_feti ! feti librairy627 !! * Local declarations628 INTEGER :: ji, jj, jk, jl629 INTEGER :: iimask(jpi,jpj)630 INTEGER :: iiend, ijend631 REAL(wp) :: zres, zres2, zdemi632 !!---------------------------------------------------------------------633 634 ! Matrix computation635 ! ------------------636 637 CALL feti_creadr(malxm,malxmax,nxm,nmorse,maan,'matrice a')638 639 nnitot = nni640 641 CALL mpp_sum( nnitot, 1, numit0ete )642 CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae')643 644 ! initialisation of the local barotropic matrix645 ! local boundary conditions on the halo646 647 CALL lbc_lnk( gcp(:,:,1), 'F', 1)648 CALL lbc_lnk( gcp(:,:,2), 'F', 1)649 CALL lbc_lnk( gcp(:,:,3), 'F', 1)650 CALL lbc_lnk( gcp(:,:,4), 'F', 1)651 CALL lbc_lnk( gcdmat , 'T', 1)652 653 ! Neumann conditions654 ! initialisation of the integer Neumann Mask655 656 CALL feti_iclr(jpi*jpj,iimask)657 DO jj = 1, jpj658 DO ji = 1, jpi659 iimask(ji,jj) = INT( gcdprc(ji,jj) )660 END DO661 END DO662 663 ! regularization of the local matrix664 665 DO jj = 1, jpj666 DO ji = 1, jpi667 gcdmat(ji,jj) = gcdmat(ji,jj) * gcdprc(ji,jj) + 1. - gcdprc(ji,jj)668 END DO669 END DO670 671 DO jk = 1, 4672 DO jj = 1, jpj673 DO ji = 1, jpi674 gcp(ji,jj,jk) = gcp(ji,jj,jk) * gcdprc(ji,jj)675 END DO676 END DO677 END DO678 679 ! implementation of the west, east, north & south Neumann conditions680 681 zdemi = 0.5682 683 ! pb with periodicity conditions : iiend, ijend684 685 ijend = nlej686 iiend = nlei687 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) iiend = nlci - jpreci688 689 IF( nbondi == 2 .AND. (nperio /= 1 .OR. nperio /= 4 .OR. nperio == 6) ) THEN690 691 ! with the periodicity : no east/west interface if nbondi = 2692 ! and nperio != 1693 694 ELSE695 ! west696 IF( nbondi /= -1 ) THEN697 DO jj = 1, jpj698 IF( iimask(1,jj) /= 0 ) THEN699 gcp(1,jj,2) = 0.e0700 gcp(1,jj,1) = zdemi * gcp(1,jj,1)701 gcp(1,jj,4) = zdemi * gcp(1,jj,4)702 ENDIF703 END DO704 DO jj = 1, jpj705 IF( iimask(1,jj) /= 0 ) THEN706 gcdmat(1,jj) = - ( gcp(1,jj,1) + gcp(1,jj,2) + gcp(1,jj,3) + gcp(1,jj,4) )707 ENDIF708 END DO709 ENDIF710 ! east711 IF( nbondi /= 1 ) THEN712 DO jj = 1, jpj713 IF( iimask(iiend,jj) /= 0 ) THEN714 gcp(iiend,jj,3) = 0.e0715 gcp(iiend,jj,1) = zdemi * gcp(iiend,jj,1)716 gcp(iiend,jj,4) = zdemi * gcp(iiend,jj,4)717 ENDIF718 END DO719 DO jj = 1, jpj720 IF( iimask(iiend,jj) /= 0 ) THEN721 gcdmat(iiend,jj) = - ( gcp(iiend,jj,1) + gcp(iiend,jj,2) &722 + gcp(iiend,jj,3) + gcp(iiend,jj,4) )723 ENDIF724 END DO725 ENDIF726 ENDIF727 728 ! south729 IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN730 DO ji = 1, jpi731 IF( iimask(ji,1) /= 0 ) THEN732 gcp(ji,1,1) = 0.e0733 gcp(ji,1,2) = zdemi * gcp(ji,1,2)734 gcp(ji,1,3) = zdemi * gcp(ji,1,3)735 ENDIF736 END DO737 DO ji = 1, jpi738 IF( iimask(ji,1) /= 0 ) THEN739 gcdmat(ji,1) = - ( gcp(ji,1,1) + gcp(ji,1,2) + gcp(ji,1,3) + gcp(ji,1,4) )740 ENDIF741 END DO742 ENDIF743 744 ! north745 IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN746 DO ji = 1, jpi747 IF( iimask(ji,ijend) /= 0 ) THEN748 gcp(ji,ijend,4) = 0.e0749 gcp(ji,ijend,2) = zdemi * gcp(ji,ijend,2)750 gcp(ji,ijend,3) = zdemi * gcp(ji,ijend,3)751 ENDIF752 END DO753 DO ji = 1, jpi754 IF( iimask(ji,ijend) /= 0 ) THEN755 gcdmat(ji,ijend) = - ( gcp(ji,ijend,1) + gcp(ji,ijend,2) &756 + gcp(ji,ijend,3) + gcp(ji,ijend,4) )757 ENDIF758 END DO759 ENDIF760 761 ! matrix terms are saved in FETI solver arrays762 CALL feti_vmov(noeuds,gcp(1,1,1),wfeti(maan))763 CALL feti_vmov(noeuds,gcp(1,1,2),wfeti(maan+noeuds))764 CALL feti_vmov(noeuds,gcdmat,wfeti(maan+2*noeuds))765 CALL feti_vmov(noeuds,gcp(1,1,3),wfeti(maan+3*noeuds))766 CALL feti_vmov(noeuds,gcp(1,1,4),wfeti(maan+4*noeuds))767 768 ! construction of Dirichlet liberty degrees array769 CALL feti_subdir(nifmat,njfmat,noeuds,ndir,iimask)770 CALL feti_creadr(malim,malimax,nim,ndir,malisdir,'lisdir')771 CALL feti_listdir(jpi,jpj,iimask,ndir,mfet(malisdir))772 773 ! stop onto matrix term for Dirichlet conditions774 CALL feti_blomat(nifmat+1,njfmat+1,wfeti(maan),ndir,mfet(malisdir))775 776 ! reservation of factorized diagonal blocs and temporary array for777 ! factorization778 npblo = (njfmat+1) * (nifmat+1) * (nifmat+1)779 ndimax = nifmat+1780 781 CALL feti_creadr(malxm,malxmax,nxm,npblo,mablo,'blo')782 CALL feti_creadr(malxm,malxmax,nxm,noeuds,madia,'dia')783 CALL feti_creadr(malxm,malxmax,nxm,noeuds,mav,'v')784 CALL feti_creadr(malxm,malxmax,nxm,ndimax*ndimax,mautil,'util')785 786 ! stoping the rigid modes787 788 ! the number of rigid modes =< Max [dim(Ker(Ep))]789 ! p=1,Np790 791 CALL feti_creadr(malim,malimax,nim,ndkerep,malisblo,'lisblo')792 793 ! Matrix factorization794 795 CALL feti_front(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo, &796 wfeti(mablo),wfeti(madia), &797 wfeti(mautil),wfeti(mav),ndlblo,mfet(malisblo),ndkerep)798 CALL feti_prext(noeuds,wfeti(madia))799 800 ! virtual dealloc => we have to see for a light f90 version801 ! the super structure is removed to clean the coarse grid802 ! solver structure803 804 malxm = madia805 CALL feti_vclr(noeuds,wfeti(madia))806 CALL feti_vclr(noeuds,wfeti(mav))807 CALL feti_vclr(ndimax*ndimax,wfeti(mautil))808 809 ! ndlblo is the dimension of the local nullspace .=<. the size of the810 ! memory of the superstructure associated to the nullspace : ndkerep811 ! ndkerep is introduced to avoid messages "out of bounds" when memory812 ! is checked813 814 ! copy matrix for Dirichlet condition815 816 CALL feti_creadr(malxm,malxmax,nxm,noeuds,miax,'x')817 CALL feti_creadr(malxm,malxmax,nxm,noeuds,may,'y')818 CALL feti_creadr(malxm,malxmax,nxm,noeuds,maz,'z')819 820 ! stoping the rigid modes821 822 ! ndlblo is the dimension of the local nullspace .=<. the size of the823 ! memory of the superstructure associated to the nullspace : ndkerep824 ! ndkerep is introduced to avoid messages "out of bounds" when memory825 ! is checked826 827 CALL feti_creadr(malxm,malxmax,nxm,ndkerep*noeuds,mansp,'nsp')828 CALL feti_blomat1(nifmat+1,njfmat+1,wfeti(maan),ndlblo, &829 mfet(malisblo),wfeti(mansp))830 831 ! computation of operator kernel832 833 CALL feti_nullsp(noeuds,nifmat+1,njfmat+1,npblo,wfeti(mablo), &834 wfeti(maan),ndlblo,mfet(malisblo),wfeti(mansp), &835 wfeti(maz))836 837 ! test of the factorisation onto each sub domain838 839 CALL feti_init(noeuds,wfeti(may))840 CALL feti_blodir(noeuds,wfeti(may),ndir,mfet(malisdir))841 CALL feti_blodir(noeuds,wfeti(may),ndlblo,mfet(malisblo))842 CALL feti_vclr(noeuds,wfeti(miax))843 CALL feti_resloc(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo, &844 wfeti(mablo),wfeti(may),wfeti(miax),wfeti(maz))845 CALL feti_proax(noeuds,nifmat+1,njfmat+1,wfeti(maan),wfeti(miax), &846 wfeti(maz))847 CALL feti_blodir(noeuds,wfeti(maz),ndlblo,mfet(malisblo))848 CALL feti_vsub(noeuds,wfeti(may),wfeti(maz),wfeti(maz))849 850 zres2 = 0.e0851 DO jl = 1, noeuds852 zres2 = zres2 + wfeti(may+jl-1) * wfeti(may+jl-1)853 END DO854 CALL mpp_sum(zres2,1,zres)855 856 res2 = 0.e0857 DO jl = 1, noeuds858 res2 = res2 + wfeti(maz+jl-1) * wfeti(maz+jl-1)859 END DO860 res2 = res2 / zres2861 CALL mpp_sum(res2,1,zres)862 863 res2 = SQRT(res2)864 IF(lwp) WRITE(numout,*) 'global residu : sqrt((Ax-b,Ax-b)/(b.b)) =', res2865 866 IF( res2 > (eps/100.) ) THEN867 IF(lwp) WRITE (numout,*) 'eps is :',eps868 IF(lwp) WRITE (numout,*) 'factorized matrix precision :',res2869 STOP870 ENDIF871 872 END SUBROUTINE fetmat873 874 875 SUBROUTINE fetsch876 !!---------------------------------------------------------------------877 !! *** ROUTINE fetsch ***878 !!879 !! ** Purpose :880 !! Construction of the matrix of the barotropic stream function881 !! system.882 !! Finite Elements Tearing & Interconnecting (FETI) approach883 !! Data framework for the Schur Dual solve884 !!885 !! ** Method :886 !!887 !! References :888 !! Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :889 !! A domain decomposition solver to compute the barotropic890 !! component of an OGCM in the parallel processing field.891 !! Ocean Modelling, issue 105, december 94.892 !!893 !! History :894 !! ! 98-02 (M. Guyon) Original code895 !! 8.5 ! 02-09 (G. Madec) F90: Free form and module896 !!----------------------------------------------------------------------897 !! * Modules used898 USE lib_feti ! feti librairy899 !! * Local declarations900 !!---------------------------------------------------------------------901 902 ! computing weights for the conform construction903 904 CALL feti_creadr(malxm,malxmax,nxm,noeuds,mapoids ,'poids' )905 CALL feti_creadr(malxm,malxmax,nxm,nnic ,mabufin ,'bufin' )906 CALL feti_creadr(malxm,malxmax,nxm,nnic ,mabufout,'bufout')907 908 !! CALL feti_poids(ninterfc,mfet(mandvoisc),mfet(maplistin),nnic, &909 !! mfet(malistin),narea,noeuds,wfeti(mapoids),wfeti(mabufin), &910 !! wfeti(mabufout) )911 CALL feti_poids(ninterfc, nnic, &912 mfet(malistin), noeuds,wfeti(mapoids) )913 914 915 ! Schur dual arrays916 917 CALL feti_creadr(malxm,malxmax,nxm,noeuds,mabitw,'bitw')918 CALL feti_creadr(malxm,malxmax,nxm,noeuds,mautilu,'utilu')919 CALL feti_creadr(malxm,malxmax,nxm,nni,malambda,'lambda')920 CALL feti_creadr(malxm,malxmax,nxm,nni,mag,'g')921 CALL feti_creadr(malxm,malxmax,nxm,nni,mapg,'pg')922 CALL feti_creadr(malxm,malxmax,nxm,nni,mamg,'mg')923 CALL feti_creadr(malxm,malxmax,nxm,nni,maw,'w')924 CALL feti_creadr(malxm,malxmax,nxm,nni,madw,'dw')925 926 ! coarse grid solver dimension and arrays927 928 nitmaxete = ndlblo929 CALL mpp_sum(nitmaxete,1,numit0ete)930 931 nitmaxete = nitmaxete + 1932 CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maxnul,'xnul')933 CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maynul,'ynul')934 CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteg,'eteg')935 CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteag,'eteag')936 CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maeted,'eted')937 CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maetead,'etead')938 CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maeteadd,'eteadd')939 CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maetegamm,'etegamm')940 CALL feti_creadr(malxm,malxmax,nxm,nni,maetew,'etew')941 CALL feti_creadr(malxm,malxmax,nxm,noeuds,maetev,'etev')942 943 ! construction of semi interface arrays944 945 CALL feti_creadr(malim,malimax,nim,ninterf+1,maplistih,'plistih')946 !! CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),nni, &947 !! mfet(maplistih),nnih,narea)948 CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin), &949 mfet(maplistih),nnih )950 951 CALL feti_creadr(malxm,malxmax,nxm,nnih,magh,'gh')952 953 ! Schur Dual Method954 955 nmaxd = nnitot / 2956 957 ! computation of the remain array for descent directions958 959 nmaxd = min0(nmaxd,(nxm-nitmaxete-malxm)/(2*nnih+3))960 CALL mpp_min(nmaxd,1,numit0ete)961 962 nitmax = nnitot/2963 epsilo = eps964 ntest = 0965 966 ! Krylov space construction967 968 CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,mawj,'wj')969 CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,madwj,'dwj')970 CALL feti_creadr(malxm,malxmax,nxm,nmaxd,madwwj,'dwwj')971 CALL feti_creadr(malxm,malxmax,nxm,nmaxd,magamm,'gamm')972 CALL feti_creadr(malxm,malxmax,nxm,max0(nmaxd,nitmaxete),mawork,'work')973 mjj0 = 0974 numit0ete = 0975 976 END SUBROUTINE fetsch977 978 #else979 SUBROUTINE fetstr ! Empty routine980 END SUBROUTINE fetstr981 SUBROUTINE fetmat ! Empty routine982 END SUBROUTINE fetmat983 SUBROUTINE fetsch ! Empty routine984 END SUBROUTINE fetsch985 #endif986 399 987 400 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
