Changeset 2587 for branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/SOL/solsor_tam.F90
- Timestamp:
- 2011-02-15T12:58:59+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/SOL/solsor_tam.F90
r1885 r2587 19 19 & jpiglo 20 20 USE in_out_manager, ONLY: & ! I/O manager 21 & nit000 21 & nit000, lwp 22 22 USE sol_oce , ONLY: & ! solver variables 23 23 & gcdmat, & … … 71 71 USE tstool_tam , ONLY: & 72 72 & prntst_adj, & 73 & stdgc 73 & stdgc, & 74 & prntst_tlm 74 75 75 76 … … 80 81 PUBLIC sol_sor_adj ! 81 82 PUBLIC sol_sor_tan ! 82 PUBLIC sol_sor_adj_tst ! called by tst.F90 83 83 PUBLIC sol_sor_adj_tst ! called by tamtst.F90 84 #if defined key_tst_tlm 85 PUBLIC sol_sor_tlm_tst ! called by tamtst.F90 86 #endif 84 87 85 88 CONTAINS … … 194 197 195 198 ! test of convergence 196 IF ( jn > nmin .AND. MOD( jn-nmin, nmod ) == 0) THEN199 IF ( (jn > nmin .AND. MOD( jn-nmin, nmod ) == 0) .OR. jn==nmax) THEN 197 200 198 201 SELECT CASE ( nsol_arp ) … … 232 235 ENDIF 233 236 ! indicator of non-convergence or explosion 237 IF( jn == nmax ) nitsor(istp) = jn 234 238 IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2 235 239 IF( ncut == 999 ) GOTO 999 … … 318 322 ijmppodd = MOD(nimpp+njmpp+jpr2di+jpr2dj+1,2) 319 323 ijpr2d = MAX(jpr2di,jpr2dj) 320 icount = 0321 324 322 325 ! Fixed number of iterations 323 326 istp = kt - nit000 + 1 324 327 iter = nitsor(istp) 325 328 icount = iter * 2 326 329 ! Output in gcx_ad 327 330 ! ---------------- … … 330 333 331 334 ! ! ============== 332 DO jn = 1, iter! Iterative loop335 DO jn = iter, 1, -1 ! Iterative loop 333 336 ! ! ============== 334 337 ! Guess red update … … 349 352 END DO 350 353 END DO 351 icount = icount + 1 352 354 icount = icount - 1 353 355 ! applied the lateral boundary conditions 354 356 IF( MOD(icount,ijpr2d+1) == 0 ) CALL lbc_lnk_e_adj( gcx_ad, c_solver_pt, 1.0_wp ) ! Lateral BCs … … 375 377 END DO 376 378 377 icount = icount + 1 378 379 icount = icount - 1 379 380 ! applied the lateral boundary conditions 380 381 IF( MOD(icount,ijpr2d+1) == 0 ) CALL lbc_lnk_e_adj( gcx_ad, c_solver_pt, 1.0_wp ) ! Lateral BCs … … 418 419 & jk, & 419 420 & kindic,& ! flags fo solver convergence 421 & kmod, & ! frequency of test for the SOR solver 420 422 & kt ! number of iteration 421 423 INTEGER, DIMENSION(jpi,jpj) :: & … … 465 467 466 468 kt=nit000 469 kindic=0 470 ! kmod = nmod ! store frequency of test for the SOR solver 471 ! nmod = 1 ! force frequency to one (remove adj_tst dependancy to nn_nmin) 467 472 468 469 473 DO jj = 1, jpj 470 474 DO ji = 1, jpi … … 491 495 END DO 492 496 END DO 493 497 ncut = 1 ! reinitilize the solver convergence flag 498 gcr_tl(:,:) = 0.0_wp 494 499 gcb_tl(:,:) = zgcb_tlin(:,:) 495 500 gcx_tl(:,:) = zgcx_tlin(:,:) … … 502 507 !-------------------------------------------------------------------- 503 508 504 DO jk = 1, jpk505 509 DO jj = nldj, nlej 506 510 DO ji = nldi, nlei … … 510 514 END DO 511 515 END DO 512 END DO513 516 !-------------------------------------------------------------------- 514 517 ! Compute the scalar product: ( L dx )^T W dy … … 520 523 ! Call the adjoint routine: dx^* = L^T dy^* 521 524 !-------------------------------------------------------------------- 522 525 gcb_ad(:,:) = 0.0_wp 523 526 gcx_ad(:,:) = zgcx_adin(:,:) 524 527 CALL sol_sor_adj(kt, kindic) … … 533 536 cl_name = 'sol_sor_adj ' 534 537 CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) 538 539 ! nmod = kmod ! restore initial frequency of test for the SOR solver 535 540 536 541 DEALLOCATE( & … … 547 552 548 553 END SUBROUTINE sol_sor_adj_tst 554 #if defined key_tst_tlm 555 SUBROUTINE sol_sor_tlm_tst( kumadt ) 556 !!----------------------------------------------------------------------- 557 !! 558 !! *** ROUTINE example_adj_tst *** 559 !! 560 !! ** Purpose : Test the tangent routine. 561 !! 562 !! ** Method : Verify the tangent with Taylor expansion 563 !! 564 !! M(x+hdx) = M(x) + L(hdx) + O(h^2) 565 !! 566 !! where L = tangent routine 567 !! M = direct routine 568 !! dx = input perturbation (random field) 569 !! h = ration on perturbation 570 !! 571 !! In the tangent test we verify that: 572 !! M(x+h*dx) - M(x) 573 !! g(h) = ------------------ ---> 1 as h ---> 0 574 !! L(h*dx) 575 !! and 576 !! g(h) - 1 577 !! f(h) = ---------- ---> k (costant) as h ---> 0 578 !! p 579 !! 580 !! History : 581 !! ! 10-02 (A. Vigilant) 582 !!----------------------------------------------------------------------- 583 #if defined key_tam 584 !! * Modules used 585 USE solsor ! Red-Black Successive Over-Relaxation solver 586 USE tamtrj ! writing out state trajectory 587 USE par_tlm, ONLY: & 588 & tlm_bch, & 589 & cur_loop, & 590 & h_ratio 591 USE istate_mod 592 USE wzvmod ! vertical velocity 593 USE gridrandom, ONLY: & 594 & grid_rd_sd 595 USE trj_tam 596 USE sol_oce , ONLY: & ! ocean dynamics and tracers variables 597 & gcb, gcx, ncut 598 USE oce , ONLY: & ! 599 & ua, ub, un 600 USE opatam_tst_ini, ONLY: & 601 & tlm_namrd 602 USE tamctl, ONLY: & ! Control parameters 603 & numtan, numtan_sc 604 !! * Arguments 605 INTEGER, INTENT(IN) :: & 606 & kumadt ! Output unit 549 607 608 !! * Local declarations 609 INTEGER :: & 610 & ji, & ! dummy loop indices 611 & jj, & 612 & jk, & 613 & kindic,& ! flags fo solver convergence 614 & kt ! number of iteration 615 INTEGER, DIMENSION(jpi,jpj) :: & 616 & iseed_2d ! 2D seed for the random number generator 617 REAL(KIND=wp) :: & 618 & zsp1, zsp2, zsp3, & ! scalar product 619 & zsp_gcb, zsp_gcx, & 620 & zsp, & 621 & gamma, & 622 & zgsp1, & 623 & zgsp2, & 624 & zgsp3, & 625 & zgsp4, & 626 & zgsp5, & 627 & zgsp6, & 628 & zgsp7 629 REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & 630 & zgcb_tlin , & ! Tangent input 631 & zgcx_tlin , & ! Tangent input 632 & zgcb_out , & ! Direct output 633 & zgcx_out , & ! Direct output 634 & zgcb_wop , & ! Direct output without perturbation 635 & zgcx_wop , & ! Direct output without perturbation 636 & zr ! 3D random field 637 CHARACTER(LEN=14) :: cl_name 638 CHARACTER (LEN=128) :: file_out, file_wop, file_xdx 639 CHARACTER (LEN=90) :: FMT 640 REAL(KIND=wp), DIMENSION(100):: & 641 & zscgcb,zscgcx, & 642 & zscerrgcb, zscerrgcx 643 INTEGER, DIMENSION(100):: & 644 & iiposgcb, ijposgcb, & 645 & iiposgcx, ijposgcx 646 INTEGER:: & 647 & ii, & 648 & isamp=40, & 649 & jsamp=40, & 650 & numsctlm 651 REAL(KIND=wp), DIMENSION(jpi,jpj) :: & 652 & zerrgcb, zerrgcx 653 654 ! Allocate memory 655 656 ALLOCATE( & 657 & zgcb_tlin( jpi,jpj), & 658 & zgcx_tlin( jpi,jpj), & 659 & zgcb_out ( jpi,jpj), & 660 & zgcx_out ( jpi,jpj), & 661 & zgcb_wop ( jpi,jpj), & 662 & zgcx_wop ( jpi,jpj), & 663 & zr( jpi,jpj) & 664 & ) 665 !================================================================== 666 ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and 667 ! dy = ( hdivb_tl, hdivn_tl ) 668 !================================================================== 669 670 !-------------------------------------------------------------------- 671 ! Reset the tangent and adjoint variables 672 !-------------------------------------------------------------------- 673 zgcb_tlin( :,:) = 0.0_wp 674 zgcx_tlin( :,:) = 0.0_wp 675 zgcb_out ( :,:) = 0.0_wp 676 zgcx_out ( :,:) = 0.0_wp 677 zgcb_wop ( :,:) = 0.0_wp 678 zgcx_wop ( :,:) = 0.0_wp 679 zr( :,:) = 0.0_wp 680 681 !-------------------------------------------------------------------- 682 ! Initialize the tangent input with random noise: dx 683 !-------------------------------------------------------------------- 684 685 !-------------------------------------------------------------------- 686 ! Output filename Xn=F(X0) 687 !-------------------------------------------------------------------- 688 CALL tlm_namrd 689 gamma = h_ratio 690 file_wop='trj_wop_solsor' 691 file_xdx='trj_xdx_solsor' 692 !-------------------------------------------------------------------- 693 ! Initialize the tangent input with random noise: dx 694 !-------------------------------------------------------------------- 695 IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN 696 CALL grid_rd_sd( 596035, zr, c_solver_pt, 0.0_wp, stdgc) 697 DO jj = nldj, nlej 698 DO ji = nldi, nlei 699 zgcb_tlin(ji,jj) = zr(ji,jj) 700 END DO 701 END DO 702 CALL grid_rd_sd( 264792, zr, c_solver_pt, 0.0_wp, stdgc) 703 DO jj = nldj, nlej 704 DO ji = nldi, nlei 705 zgcx_tlin(ji,jj) = zr(ji,jj) 706 END DO 707 END DO 708 ENDIF 709 710 !-------------------------------------------------------------------- 711 ! Complete Init for Direct 712 !------------------------------------------------------------------- 713 CALL istate_p 714 715 ! *** initialize the reference trajectory 716 ! ------------ 717 718 ! gcx (:,:) = ( ua(:,:,1) + ub(:,:,1) ) / 10.0_wp 719 ! gcb (:,:) = ( ua(:,:,3) + ub(:,:,3) ) / 10.0_wp 720 CALL trj_rea( nit000-1, 1 ) 721 CALL trj_rea( nit000, 1 ) 722 gcx (:,:) = un(:,:,1) / 10.0_wp 723 gcb (:,:) = un(:,:,3) / 10.0_wp 724 725 IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN 726 zgcb_tlin(:,:) = gamma * zgcb_tlin(:,:) 727 gcb(:,:) = gcb(:,:) + zgcb_tlin(:,:) 728 729 zgcx_tlin(:,:) = gamma * zgcx_tlin(:,:) 730 gcx(:,:) = gcx(:,:) + zgcx_tlin(:,:) 731 ENDIF 732 733 !-------------------------------------------------------------------- 734 ! Compute the direct model F(X0,t=n) = Xn 735 !-------------------------------------------------------------------- 736 kindic=0 737 ncut=1 738 IF ( tlm_bch /= 2 ) CALL sol_sor(kindic) 739 IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop) 740 IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx) 741 !-------------------------------------------------------------------- 742 ! Compute the Tangent 743 !-------------------------------------------------------------------- 744 IF ( tlm_bch == 2 ) THEN 745 !-------------------------------------------------------------------- 746 ! Initialize the tangent variables: dy^* = W dy 747 !-------------------------------------------------------------------- 748 gcr_tl(:,:) = 0.0_wp 749 gcb_tl (:,:) = zgcb_tlin (:,:) 750 gcx_tl (:,:) = zgcx_tlin (:,:) 751 752 !----------------------------------------------------------------------- 753 ! Initialization of the dynamics and tracer fields for the tangent 754 !----------------------------------------------------------------------- 755 ncut=1 !reset indicator of solver convergence 756 CALL sol_sor_tan(nit000, kindic) 757 !-------------------------------------------------------------------- 758 ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) 759 !-------------------------------------------------------------------- 760 761 zsp_gcx = DOT_PRODUCT( gcx_tl, gcx_tl ) 762 zsp2 = zsp_gcx 763 !-------------------------------------------------------------------- 764 ! Storing data 765 !-------------------------------------------------------------------- 766 CALL trj_rd_spl(file_wop) 767 zgcx_wop (:,:) = gcx (:,:) 768 CALL trj_rd_spl(file_xdx) 769 zgcx_out (:,:) = gcx (:,:) 770 !-------------------------------------------------------------------- 771 ! Compute the Linearization Error 772 ! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn) 773 ! and 774 ! Compute the Linearization Error 775 ! En = Nn -TL(gamma.dX0, t0,tn) 776 !-------------------------------------------------------------------- 777 ! Warning: Here we re-use local variables z()_out and z()_wop 778 ii=0 779 DO jj = 1, jpj 780 DO ji = 1, jpi 781 zgcx_out (ji,jj) = zgcx_out (ji,jj) - zgcx_wop (ji,jj) 782 zgcx_wop (ji,jj) = zgcx_out (ji,jj) - gcx_tl (ji,jj) 783 IF ( gcx_tl(ji,jj) .NE. 0.0_wp ) zerrgcx(ji,jj) = zgcx_out(ji,jj)/gcx_tl(ji,jj) 784 IF( (MOD(ji, isamp) .EQ. 0) .AND. & 785 & (MOD(jj, jsamp) .EQ. 0) ) THEN 786 ii = ii+1 787 iiposgcx(ii) = ji 788 ijposgcx(ii) = jj 789 IF ( INT(tmask(ji,jj,1)) .NE. 0) THEN 790 zscgcx (ii) = zgcx_wop(ji,jj) 791 zscerrgcx (ii) = ( zerrgcx(ji,jj) - 1.0_wp ) / gamma 792 ENDIF 793 ENDIF 794 END DO 795 END DO 796 797 zsp_gcx = DOT_PRODUCT( zgcx_out, zgcx_out ) 798 799 zsp1 = zsp_gcx 800 801 zsp_gcx = DOT_PRODUCT( zgcx_wop, zgcx_wop ) 802 803 zsp3 = zsp_gcx 804 !-------------------------------------------------------------------- 805 ! Print the linearization error En - norme 2 806 !-------------------------------------------------------------------- 807 ! 14 char:'12345678901234' 808 cl_name = 'sol_sor: En ' 809 zsp = SQRT(zsp3) 810 zgsp5 = zsp 811 CALL prntst_tlm( cl_name, kumadt, zsp, h_ratio ) 812 813 !-------------------------------------------------------------------- 814 ! Compute TLM norm2 815 !-------------------------------------------------------------------- 816 zsp = SQRT(zsp2) 817 818 zgsp4 = zsp 819 cl_name = 'sol_sor: Ln2 ' 820 CALL prntst_tlm( cl_name, kumadt, zsp, h_ratio ) 821 822 !-------------------------------------------------------------------- 823 ! Print the linearization error Nn - norme 2 824 !-------------------------------------------------------------------- 825 zsp = SQRT(zsp1) 826 827 cl_name = 'solsor:Mhdx-Mx' 828 CALL prntst_tlm( cl_name, kumadt, zsp, h_ratio ) 829 830 zgsp3 = SQRT( zsp3/zsp2 ) 831 zgsp7 = zgsp3/gamma 832 zgsp1 = zsp 833 zgsp2 = zgsp1 / zgsp4 834 zgsp6 = (zgsp2 - 1.0_wp)/gamma 835 836 FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)" 837 WRITE(numtan,FMT) 'solsor ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 838 !-------------------------------------------------------------------- 839 ! Unitary calculus 840 !-------------------------------------------------------------------- 841 FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)" 842 cl_name ='sol_sor ' 843 IF(lwp) THEN 844 DO ii=1, 100, 1 845 IF ( zscgcx(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgcx ', & 846 & cur_loop, h_ratio, ii, iiposgcx(ii), ijposgcx(ii), zscgcx(ii) 847 ENDDO 848 DO ii=1, 100, 1 849 IF ( zscerrgcx(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrgcx ', & 850 & cur_loop, h_ratio, ii, iiposgcx(ii), ijposgcx(ii), zscerrgcx(ii) 851 ENDDO 852 ! write separator 853 WRITE(numtan_sc,"(A4)") '====' 854 ENDIF 855 856 ENDIF 857 858 DEALLOCATE( & 859 & zgcb_tlin, & 860 & zgcx_tlin, & 861 & zgcb_out , & 862 & zgcx_out , & 863 & zgcb_wop , & 864 & zgcx_wop , & 865 & zr & 866 & ) 867 #else 868 !! * Arguments 869 INTEGER, INTENT(IN) :: & 870 & kumadt ! Output unit 871 ! dummy routine 872 #endif 873 END SUBROUTINE sol_sor_tlm_tst 874 #endif 550 875 END MODULE solsor_tam
Note: See TracChangeset
for help on using the changeset viewer.
