Changeset 3501 for branches/2012
- Timestamp:
- 2012-10-15T17:56:16+02:00 (11 years ago)
- Location:
- branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
r3452 r3501 52 52 53 53 !! * Routine accessibility 54 PUBLIC dia_dct ! routine called by step.F90 55 PUBLIC dia_dct_init! routine called by opa.F90 54 PUBLIC dia_dct ! routine called by step.F90 55 PUBLIC dia_dct_init ! routine called by opa.F90 56 PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90 56 57 PRIVATE readsec 57 58 PRIVATE removepoints … … 73 74 INTEGER, PARAMETER :: nb_point_max = 2000 74 75 INTEGER, PARAMETER :: nb_type_class = 14 76 INTEGER, PARAMETER :: nb_3d_vars = 5 77 INTEGER, PARAMETER :: nb_2d_vars = 2 75 78 INTEGER :: nb_sec 76 79 … … 106 109 TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 107 110 108 111 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d 112 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d 113 109 114 CONTAINS 115 116 117 INTEGER FUNCTION diadct_alloc() 118 !!---------------------------------------------------------------------- 119 !! *** FUNCTION diadct_alloc *** 120 !!---------------------------------------------------------------------- 121 INTEGER :: ierr(2) 122 !!---------------------------------------------------------------------- 123 124 ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 125 ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(2) ) 126 127 diadct_alloc = MAXVAL( ierr ) 128 IF( diadct_alloc /= 0 ) CALL ctl_warn('diadct_alloc: failed to allocate arrays') 129 130 END FUNCTION diadct_alloc 110 131 111 132 SUBROUTINE dia_dct_init … … 113 134 !! *** ROUTINE diadct *** 114 135 !! 115 !! ** Purpose: Read the namelist paramet res136 !! ** Purpose: Read the namelist parameters 116 137 !! Open output files 117 138 !! … … 154 175 ENDIF 155 176 177 ! Initialise arrays to zero 178 transports_3d(:,:,:,:)=0.0 179 transports_2d(:,:,:) =0.0 180 156 181 IF( nn_timing == 1 ) CALL timing_stop('dia_dct_init') 157 182 ! … … 163 188 !! *** ROUTINE diadct *** 164 189 !! 165 !! ** Purpose: Compute sections tranport and write it in numdct file 190 !! Purpose :: Compute section transports and write it in numdct files 191 !! 192 !! Method :: All arrays initialised to zero in dct_init 193 !! Each nn_dct time step call subroutine 'transports' for 194 !! each section to sum the transports over each grid cell. 195 !! Each nn_dctwri time step: 196 !! Divide the arrays by the number of summations to gain 197 !! an average value 198 !! Call dia_dct_sum to sum relevant grid boxes to obtain 199 !! totals for each class (density, depth, temp or sal) 200 !! Call dia_dct_wri to write the transports into file 201 !! Reinitialise all relevant arrays to zero 166 202 !!--------------------------------------------------------------------- 167 203 !! * Arguments … … 170 206 !! * Local variables 171 207 INTEGER :: jsec, &! loop on sections 172 iost, &! error for opening fileout173 208 itotal ! nb_sec_max*nb_type_class*nb_class_max 174 209 LOGICAL :: lldebug =.FALSE. ! debug a section 175 CHARACTER(len=160) :: clfileout ! fileout name176 177 210 178 211 INTEGER , DIMENSION(1) :: ish ! tmp array for mpp_sum … … 208 241 209 242 !Compute transport through section 210 CALL transport(secs(jsec),lldebug )243 CALL transport(secs(jsec),lldebug,jsec) 211 244 212 245 ENDDO … … 214 247 IF( MOD(kt,nn_dctwri)==0 )THEN 215 248 216 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: write at kt = ",kt249 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average transports and write at kt = ",kt 217 250 251 !! divide arrays by nn_dctwri/nn_dct to obtain average 252 transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 253 transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct) 254 255 ! Sum over each class 256 DO jsec=1,nb_sec 257 CALL dia_dct_sum(secs(jsec),jsec) 258 ENDDO 259 218 260 !Sum on all procs 219 261 IF( lk_mpp )THEN … … 233 275 234 276 !nullify transports values after writing 277 transports_3d(:,jsec,:,:)=0. 278 transports_2d(:,jsec,: )=0. 235 279 secs(jsec)%transport(:,:)=0. 236 280 … … 535 579 END SUBROUTINE removepoints 536 580 537 SUBROUTINE transport(sec,ld_debug )581 SUBROUTINE transport(sec,ld_debug,jsec) 538 582 !!------------------------------------------------------------------------------------------- 539 583 !! *** ROUTINE transport *** 540 584 !! 541 !! ** Purpose : Compute the transport through a section 542 !! 543 !! ** Method :Transport through a given section is equal to the sum of transports 544 !! computed on each proc. 545 !! On each proc,transport is equal to the sum of transport computed through 546 !! segments linking each point of sec%listPoint with the next one. 547 !! 548 !! !BE carefull : 549 !! one section is a sum of segments 550 !! one segment is defined by 2 consectuve points in sec%listPoint 551 !! all points of sec%listPoint are positioned on the F-point of the cell. 585 !! Purpose :: Compute the transport for each point in a section 552 586 !! 553 !! There are several loops: 554 !! loop on the density/temperature/salinity/level classes 555 !! loop on the segment between 2 nodes 556 !! loop on the level jk 557 !! test on the density/temperature/salinity/level 558 !! 559 !! ** Output: sec%transport: volume/mass/ice/heat/salt transport in the 2 directions 560 !! 587 !! Method :: Loop over each segment, and each vertical level and add the transport 588 !! Be aware : 589 !! One section is a sum of segments 590 !! One segment is defined by 2 consecutive points in sec%listPoint 591 !! All points of sec%listPoint are positioned on the F-point of the cell 592 !! 593 !! There are two loops: 594 !! loop on the segment between 2 nodes 595 !! loop on the level jk !! 596 !! 597 !! Output :: Arrays containing the volume,density,salinity,temperature etc 598 !! transports for each point in a section, summed over each nn_dct. 561 599 !! 562 600 !!------------------------------------------------------------------------------------------- … … 564 602 TYPE(SECTION),INTENT(INOUT) :: sec 565 603 LOGICAL ,INTENT(IN) :: ld_debug 604 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 566 605 567 606 !! * Local variables 568 INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes 569 isgnu , isgnv ! 570 INTEGER :: ii, ij ! local integer 571 REAL(wp):: zumid , zvmid ,&!U/V velocity on a cell segment 572 zumid_ice , zvmid_ice ,&!U/V ice velocity 573 zTnorm ,&!transport of velocity through one cell's sides 574 ztransp1 , ztransp2 ,&!total transport in directions 1 and 2 575 ztemp1 , ztemp2 ,&!temperature transport " 576 zrhoi1 , zrhoi2 ,&!mass transport " 577 zrhop1 , zrhop2 ,&!mass transport " 578 zsal1 , zsal2 ,&!salinity transport " 579 zice_vol_pos , zice_vol_neg ,&!volume ice transport " 580 zice_surf_pos, zice_surf_neg !surface ice transport " 581 REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 607 INTEGER :: jk, jseg, jclass, &!loop on level/segment/classes 608 isgnu, isgnv ! 609 REAL(wp) :: zumid, zvmid, &!U/V velocity on a cell segment 610 zumid_ice, zvmid_ice, &!U/V ice velocity 611 zTnorm !transport of velocity through one cell's sides 612 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point 582 613 583 614 TYPE(POINT_SECTION) :: k 584 REAL(wp), POINTER, DIMENSION(:,:):: zsum ! 2D work array585 615 !!-------------------------------------------------------- 586 CALL wrk_alloc( nb_type_class , nb_class_max , zsum )587 616 588 617 IF( ld_debug )WRITE(numout,*)' Compute transport' 589 590 !----------------!591 ! INITIALIZATION !592 !----------------!593 zsum = 0._wp594 zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp595 zice_vol_pos = 0._wp ; zice_vol_neg = 0._wp596 618 597 619 !---------------------------! … … 670 692 END SELECT 671 693 672 !------------------------------- 673 ! LOOP ON THE DENSITY CLASSES | 674 !------------------------------- 675 !The computation is made for each density class 676 DO jclass=1,MAX(1,sec%nb_class-1) 677 678 ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp 679 ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp 680 681 !---------------------------| 682 ! LOOP ON THE LEVEL | 683 !---------------------------| 684 !Sum of the transport on the vertical 685 DO jk=1,jpk 686 687 688 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 689 SELECT CASE( sec%direction(jseg) ) 690 CASE(0,1) 691 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 692 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 693 zrhop = interp(k%I,k%J,jk,'V',rhop) 694 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 695 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 696 CASE(2,3) 697 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 698 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 699 zrhop = interp(k%I,k%J,jk,'U',rhop) 700 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 701 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 702 END SELECT 703 704 zfsdep= gdept(k%I,k%J,jk) 705 706 !----------------------------------------------! 707 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 708 !----------------------------------------------! 709 710 IF ( ( ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 711 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 712 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 713 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 714 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 715 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 716 ((( zsn .GT. sec%zsal(jclass)) .AND. & 717 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 718 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 719 ((( ztn .GE. sec%ztem(jclass)) .AND. & 720 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 721 ( sec%ztem(jclass) .EQ.99.)) .AND. & 722 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 723 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 724 ( sec%zlay(jclass) .EQ. 99. )))) THEN 725 726 727 !compute velocity with the correct direction 728 SELECT CASE( sec%direction(jseg) ) 729 CASE(0,1) 730 zumid=0. 731 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 732 CASE(2,3) 733 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 734 zvmid=0. 735 END SELECT 736 737 !velocity* cell's length * cell's thickness 738 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 739 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 694 !---------------------------| 695 ! LOOP ON THE LEVEL | 696 !---------------------------| 697 !Sum of the transport on the vertical 698 DO jk=1,mbathy(k%I,k%J) 699 700 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 701 SELECT CASE( sec%direction(jseg) ) 702 CASE(0,1) 703 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 704 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 705 zrhop = interp(k%I,k%J,jk,'V',rhop) 706 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 707 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 708 CASE(2,3) 709 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 710 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 711 zrhop = interp(k%I,k%J,jk,'U',rhop) 712 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 713 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 714 END SELECT 715 716 zfsdep= gdept(k%I,k%J,jk) 717 718 !compute velocity with the correct direction 719 SELECT CASE( sec%direction(jseg) ) 720 CASE(0,1) 721 zumid=0. 722 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 723 CASE(2,3) 724 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 725 zvmid=0. 726 END SELECT 727 728 !zTnorm=transport through one cell; 729 !velocity* cell's length * cell's thickness 730 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 731 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 740 732 741 733 #if ! defined key_vvl 742 !add transport due to free surface743 IF( jk==1 )THEN744 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &745 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)746 ENDIF734 !add transport due to free surface 735 IF( jk==1 )THEN 736 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 737 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 738 ENDIF 747 739 #endif 748 !COMPUTE TRANSPORT 749 !zTnorm=transport through one cell for one class 750 !ztransp1 or ztransp2=transport through one cell i 751 ! for one class for one direction 752 IF( zTnorm .GE. 0 )THEN 753 754 ztransp1=zTnorm+ztransp1 755 756 IF ( sec%llstrpond ) THEN 757 ztemp1 = ztemp1 + zTnorm * ztn 758 zsal1 = zsal1 + zTnorm * zsn 759 zrhoi1 = zrhoi1 + zTnorm * zrhoi 760 zrhop1 = zrhop1 + zTnorm * zrhop 761 ENDIF 762 763 ELSE 764 765 ztransp2=(zTnorm)+ztransp2 766 767 IF ( sec%llstrpond ) THEN 768 ztemp2 = ztemp2 + zTnorm * ztn 769 zsal2 = zsal2 + zTnorm * zsn 770 zrhoi2 = zrhoi2 + zTnorm * zrhoi 771 zrhop2 = zrhop2 + zTnorm * zrhop 772 ENDIF 773 ENDIF 774 775 776 ENDIF ! end of density test 777 ENDDO!end of loop on the level 778 779 !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE DIRECTIONS 780 !--------------------------------------------------- 781 zsum(1,jclass) = zsum(1,jclass)+ztransp1 782 zsum(2,jclass) = zsum(2,jclass)+ztransp2 783 IF( sec%llstrpond )THEN 784 zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1 785 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2 786 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1 787 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2 788 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1 789 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2 790 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1 791 zsum(10,jclass) = zsum(10,jclass)+zsal2 740 !COMPUTE TRANSPORT 741 742 transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 743 744 IF ( sec%llstrpond ) THEN 745 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * zrhoi 746 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zrhop 747 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk) + zTnorm * ztn 748 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk) + zTnorm * zsn 792 749 ENDIF 793 750 794 ENDDO !end of loop on the density classes751 ENDDO !end of loop on the level 795 752 796 753 #if defined key_lim2 || defined key_lim3 … … 816 773 zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 817 774 818 IF( zTnorm .GE. 0)THEN 819 zice_vol_pos = (zTnorm)* & 820 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 821 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 822 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 823 +zice_vol_pos 824 zice_surf_pos = (zTnorm)* & 825 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 826 +zice_surf_pos 827 ELSE 828 zice_vol_neg=(zTnorm)* & 829 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 830 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 831 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 832 +zice_vol_neg 833 zice_surf_neg=(zTnorm)* & 834 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 835 +zice_surf_neg 836 ENDIF 837 838 zsum(11,1) = zsum(11,1)+zice_vol_pos 839 zsum(12,1) = zsum(12,1)+zice_vol_neg 840 zsum(13,1) = zsum(13,1)+zice_surf_pos 841 zsum(14,1) = zsum(14,1)+zice_surf_neg 775 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* & 776 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 777 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 778 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 779 +zice_vol_pos 780 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* & 781 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 782 +zice_surf_pos 842 783 843 784 ENDIF !end of ice case … … 846 787 ENDDO !end of loop on the segment 847 788 848 849 ELSE !if sec%nb_point =0850 zsum(1:2,:)=0.851 IF (sec%llstrpond) zsum(3:10,:)=0.852 zsum( 11:14,:)=0.853 789 ENDIF !end of sec%nb_point =0 case 854 855 !-------------------------------|856 !FINISH COMPUTING TRANSPORTS |857 !-------------------------------|858 DO jclass=1,MAX(1,sec%nb_class-1)859 sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6860 sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6861 IF( sec%llstrpond ) THEN862 IF( zsum(1,jclass) .NE. 0._wp ) THEN863 sec%transport( 3,jclass) = sec%transport( 3,jclass) + zsum( 3,jclass)/zsum(1,jclass)864 sec%transport( 5,jclass) = sec%transport( 5,jclass) + zsum( 5,jclass)/zsum(1,jclass)865 sec%transport( 7,jclass) = sec%transport( 7,jclass) + zsum( 7,jclass)866 sec%transport( 9,jclass) = sec%transport( 9,jclass) + zsum( 9,jclass)867 ENDIF868 IF( zsum(2,jclass) .NE. 0._wp )THEN869 sec%transport( 4,jclass) = sec%transport( 4,jclass) + zsum( 4,jclass)/zsum(2,jclass)870 sec%transport( 6,jclass) = sec%transport( 6,jclass) + zsum( 6,jclass)/zsum(2,jclass)871 sec%transport( 8,jclass) = sec%transport( 8,jclass) + zsum( 8,jclass)872 sec%transport(10,jclass) = sec%transport(10,jclass) + zsum(10,jclass)873 ENDIF874 ELSE875 sec%transport( 3,jclass) = 0._wp876 sec%transport( 4,jclass) = 0._wp877 sec%transport( 5,jclass) = 0._wp878 sec%transport( 6,jclass) = 0._wp879 sec%transport( 7,jclass) = 0._wp880 sec%transport( 8,jclass) = 0._wp881 sec%transport(10,jclass) = 0._wp882 ENDIF883 ENDDO884 885 IF( sec%ll_ice_section ) THEN886 sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6887 sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6888 sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6889 sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6890 ENDIF891 892 CALL wrk_dealloc( nb_type_class , nb_class_max , zsum )893 790 ! 894 791 END SUBROUTINE transport 792 793 SUBROUTINE dia_dct_sum(sec,jsec) 794 !!------------------------------------------------------------- 795 !! Purpose: Average the transport over nn_dctwri time steps 796 !! and sum over the density/salinity/temperature/depth classes 797 !! 798 !! Method: 799 !! Sum over relevant grid cells to obtain values 800 !! for each 801 !! There are several loops: 802 !! loop on the segment between 2 nodes 803 !! loop on the level jk 804 !! loop on the density/temperature/salinity/level classes 805 !! test on the density/temperature/salinity/level 806 !! 807 !! ** Method :Transport through a given section is equal to the sum of transports 808 !! computed on each proc. 809 !! On each proc,transport is equal to the sum of transport computed through 810 !! segments linking each point of sec%listPoint with the next one. 811 !! 812 !!------------------------------------------------------------- 813 !! * arguments 814 TYPE(SECTION),INTENT(INOUT) :: sec 815 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 816 817 TYPE(POINT_SECTION) :: k 818 INTEGER :: jk,jseg,jclass !loop on level/segment/classes 819 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 820 !!------------------------------------------------------------- 821 822 !! Sum the relevant segments to obtain values for each class 823 IF(sec%nb_point .NE. 0)THEN 824 825 !--------------------------------------! 826 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 827 !--------------------------------------! 828 DO jseg=1,MAX(sec%nb_point-1,0) 829 830 !------------------------------------------------------------------------------------------- 831 ! Select the appropriate coordinate for computing the velocity of the segment 832 ! 833 ! CASE(0) Case (2) 834 ! ------- -------- 835 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 836 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 837 ! | 838 ! | 839 ! | 840 ! Case (3) U(i,j) 841 ! -------- | 842 ! | 843 ! listPoint(jseg+1) F(i,j+1) | 844 ! | | 845 ! | | 846 ! | listPoint(jseg+1) F(i,j-1) 847 ! | 848 ! | 849 ! U(i,j+1) 850 ! | Case(1) 851 ! | ------ 852 ! | 853 ! | listPoint(jseg+1) listPoint(jseg) 854 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 855 ! listPoint(jseg) F(i,j) 856 ! 857 !------------------------------------------------------------------------------------------- 858 859 SELECT CASE( sec%direction(jseg) ) 860 CASE(0) ; k = sec%listPoint(jseg) 861 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 862 CASE(2) ; k = sec%listPoint(jseg) 863 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 864 END SELECT 865 866 !---------------------------| 867 ! LOOP ON THE LEVEL | 868 !---------------------------| 869 !Sum of the transport on the vertical 870 DO jk=1,mbathy(k%I,k%J) 871 872 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 873 SELECT CASE( sec%direction(jseg) ) 874 CASE(0,1) 875 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 876 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 877 zrhop = interp(k%I,k%J,jk,'V',rhop) 878 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 879 880 CASE(2,3) 881 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 882 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 883 zrhop = interp(k%I,k%J,jk,'U',rhop) 884 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 885 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 886 END SELECT 887 888 zfsdep= gdept(k%I,k%J,jk) 889 890 !------------------------------- 891 ! LOOP ON THE DENSITY CLASSES | 892 !------------------------------- 893 !The computation is made for each density/heat/salt/... class 894 DO jclass=1,MAX(1,sec%nb_class-1) 895 896 !----------------------------------------------! 897 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 898 !----------------------------------------------! 899 900 IF ( ( & 901 ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 902 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 903 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 904 905 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 906 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 907 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 908 909 ((( zsn .GT. sec%zsal(jclass)) .AND. & 910 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 911 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 912 913 ((( ztn .GE. sec%ztem(jclass)) .AND. & 914 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 915 ( sec%ztem(jclass) .EQ.99.)) .AND. & 916 917 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 918 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 919 ( sec%zlay(jclass) .EQ. 99. )) & 920 )) THEN 921 922 !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 923 !---------------------------------------------------------------------------- 924 IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 925 sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 926 ELSE 927 sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 928 ENDIF 929 IF( sec%llstrpond )THEN 930 931 IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN 932 933 IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 934 sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 935 ELSE 936 sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 937 ENDIF 938 939 IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 940 sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 941 ELSE 942 sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 943 ENDIF 944 945 ENDIF 946 947 IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 948 sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 949 ELSE 950 sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 951 ENDIF 952 953 IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 954 sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 955 ELSE 956 sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 957 ENDIF 958 959 ELSE 960 sec%transport( 3,jclass) = 0._wp 961 sec%transport( 4,jclass) = 0._wp 962 sec%transport( 5,jclass) = 0._wp 963 sec%transport( 6,jclass) = 0._wp 964 sec%transport( 7,jclass) = 0._wp 965 sec%transport( 8,jclass) = 0._wp 966 sec%transport( 9,jclass) = 0._wp 967 sec%transport(10,jclass) = 0._wp 968 ENDIF 969 970 ENDIF ! end of test if point is in class 971 972 ENDDO ! end of loop on the classes 973 974 ENDDO ! loop over jk 975 976 #if defined key_lim2 || defined key_lim3 977 978 !ICE CASE 979 IF( sec%ll_ice_section )THEN 980 981 IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 982 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6 983 ELSE 984 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6 985 ENDIF 986 987 IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 988 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6 989 ELSE 990 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6 991 ENDIF 992 993 ENDIF !end of ice case 994 #endif 995 996 ENDDO !end of loop on the segment 997 998 ELSE !if sec%nb_point =0 999 sec%transport(1:2,:)=0. 1000 IF (sec%llstrpond) sec%transport(3:10,:)=0. 1001 IF (sec%ll_ice_section) sec%transport( 11:14,:)=0. 1002 ENDIF !end of sec%nb_point =0 case 1003 1004 END SUBROUTINE dia_dct_sum 895 1005 896 1006 SUBROUTINE dia_dct_wri(kt,ksec,sec) … … 917 1027 918 1028 !!local declarations 919 INTEGER :: jcl ,ji! Dummy loop1029 INTEGER :: jclass ! Dummy loop 920 1030 CHARACTER(len=2) :: classe ! Classname 921 1031 REAL(wp) :: zbnd1,zbnd2 ! Class bounds 922 1032 REAL(wp) :: zslope ! section's slope coeff 923 1033 ! 924 REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace1034 REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 925 1035 !!------------------------------------------------------------- 926 CALL wrk_alloc(nb_type_class , zsumclass )927 928 zsumclass (:)=0._wp1036 CALL wrk_alloc(nb_type_class , zsumclasses ) 1037 1038 zsumclasses(:)=0._wp 929 1039 zslope = sec%slopeSection 930 1040 931 1041 932 DO jcl=1,MAX(1,sec%nb_class-1) 933 934 ! Mean computation 935 sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct) 1042 DO jclass=1,MAX(1,sec%nb_class-1) 1043 936 1044 classe = 'N ' 937 1045 zbnd1 = 0._wp 938 1046 zbnd2 = 0._wp 939 zsumclass (1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)1047 zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass) 940 1048 941 1049 942 1050 !insitu density classes transports 943 IF( ( sec%zsigi(jcl ) .NE. 99._wp ) .AND. &944 ( sec%zsigi(jcl +1) .NE. 99._wp ) )THEN1051 IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. & 1052 ( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN 945 1053 classe = 'DI ' 946 zbnd1 = sec%zsigi(jcl )947 zbnd2 = sec%zsigi(jcl +1)1054 zbnd1 = sec%zsigi(jclass) 1055 zbnd2 = sec%zsigi(jclass+1) 948 1056 ENDIF 949 1057 !potential density classes transports 950 IF( ( sec%zsigp(jcl ) .NE. 99._wp ) .AND. &951 ( sec%zsigp(jcl +1) .NE. 99._wp ) )THEN1058 IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. & 1059 ( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN 952 1060 classe = 'DP ' 953 zbnd1 = sec%zsigp(jcl )954 zbnd2 = sec%zsigp(jcl +1)1061 zbnd1 = sec%zsigp(jclass) 1062 zbnd2 = sec%zsigp(jclass+1) 955 1063 ENDIF 956 1064 !depth classes transports 957 IF( ( sec%zlay(jcl ) .NE. 99._wp ) .AND. &958 ( sec%zlay(jcl +1) .NE. 99._wp ) )THEN1065 IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. & 1066 ( sec%zlay(jclass+1) .NE. 99._wp ) )THEN 959 1067 classe = 'Z ' 960 zbnd1 = sec%zlay(jcl )961 zbnd2 = sec%zlay(jcl +1)1068 zbnd1 = sec%zlay(jclass) 1069 zbnd2 = sec%zlay(jclass+1) 962 1070 ENDIF 963 1071 !salinity classes transports 964 IF( ( sec%zsal(jcl ) .NE. 99._wp ) .AND. &965 ( sec%zsal(jcl +1) .NE. 99._wp ) )THEN1072 IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. & 1073 ( sec%zsal(jclass+1) .NE. 99._wp ) )THEN 966 1074 classe = 'S ' 967 zbnd1 = sec%zsal(jcl )968 zbnd2 = sec%zsal(jcl +1)1075 zbnd1 = sec%zsal(jclass) 1076 zbnd2 = sec%zsal(jclass+1) 969 1077 ENDIF 970 1078 !temperature classes transports 971 IF( ( sec%ztem(jcl ) .NE. 99._wp ) .AND. &972 ( sec%ztem(jcl +1) .NE. 99._wp ) ) THEN1079 IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. & 1080 ( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN 973 1081 classe = 'T ' 974 zbnd1 = sec%ztem(jcl )975 zbnd2 = sec%ztem(jcl +1)1082 zbnd1 = sec%ztem(jclass) 1083 zbnd2 = sec%ztem(jclass+1) 976 1084 ENDIF 977 1085 978 1086 !write volume transport per class 979 1087 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 980 jcl ,classe,zbnd1,zbnd2,&981 sec%transport(1,jcl ),sec%transport(2,jcl), &982 sec%transport(1,jcl )+sec%transport(2,jcl)1088 jclass,classe,zbnd1,zbnd2,& 1089 sec%transport(1,jclass),sec%transport(2,jclass), & 1090 sec%transport(1,jclass)+sec%transport(2,jclass) 983 1091 984 1092 IF( sec%llstrpond )THEN … … 986 1094 !write heat transport per class: 987 1095 WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 988 jcl ,classe,zbnd1,zbnd2,&989 sec%transport(7,jcl )*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &990 ( sec%transport(7,jcl )+sec%transport(8,jcl) )*1000._wp*rcp/1.e151096 jclass,classe,zbnd1,zbnd2,& 1097 sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, & 1098 ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15 991 1099 !write salt transport per class 992 1100 WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 993 jcl ,classe,zbnd1,zbnd2,&994 sec%transport(9,jcl )*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&995 (sec%transport(9,jcl )+sec%transport(10,jcl))*1000._wp/1.e91101 jclass,classe,zbnd1,zbnd2,& 1102 sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,& 1103 (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9 996 1104 ENDIF 997 1105 … … 1000 1108 zbnd1 = 0._wp 1001 1109 zbnd2 = 0._wp 1002 jcl =01110 jclass=0 1003 1111 1004 1112 !write total volume transport 1005 1113 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 1006 jcl ,"total",zbnd1,zbnd2,&1007 zsumclass (1),zsumclass(2),zsumclass(1)+zsumclass(2)1114 jclass,"total",zbnd1,zbnd2,& 1115 zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2) 1008 1116 1009 1117 IF( sec%llstrpond )THEN … … 1011 1119 !write total heat transport 1012 1120 WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 1013 jcl ,"total",zbnd1,zbnd2,&1014 zsumclass (7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,&1015 (zsumclass (7)+zsumclass(8) )* 1000._wp*rcp/1.e151121 jclass,"total",zbnd1,zbnd2,& 1122 zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,& 1123 (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15 1016 1124 !write total salt transport 1017 1125 WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 1018 jcl ,"total",zbnd1,zbnd2,&1019 zsumclass (9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,&1020 (zsumclass (9)+zsumclass(10))*1000._wp/1.e91126 jclass,"total",zbnd1,zbnd2,& 1127 zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,& 1128 (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9 1021 1129 ENDIF 1022 1130 … … 1025 1133 !write total ice volume transport 1026 1134 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1027 jcl ,"ice_vol",zbnd1,zbnd2,&1135 jclass,"ice_vol",zbnd1,zbnd2,& 1028 1136 sec%transport(9,1),sec%transport(10,1),& 1029 1137 sec%transport(9,1)+sec%transport(10,1) 1030 1138 !write total ice surface transport 1031 1139 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1032 jcl ,"ice_surf",zbnd1,zbnd2,&1140 jclass,"ice_surf",zbnd1,zbnd2,& 1033 1141 sec%transport(11,1),sec%transport(12,1), & 1034 1142 sec%transport(11,1)+sec%transport(12,1) … … 1038 1146 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 1039 1147 1040 CALL wrk_dealloc(nb_type_class , zsumclass )1148 CALL wrk_dealloc(nb_type_class , zsumclasses ) 1041 1149 END SUBROUTINE dia_dct_wri 1042 1150 -
branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3352 r3501 512 512 USE ldftra_oce, ONLY: ldftra_oce_alloc 513 513 USE trc_oce , ONLY: trc_oce_alloc 514 #if defined key_diadct 515 USE diadct , ONLY: diadct_alloc 516 #endif 514 517 ! 515 518 INTEGER :: ierr … … 525 528 ierr = ierr + lib_mpp_alloc (numout) ! mpp exchanges 526 529 ierr = ierr + trc_oce_alloc () ! shared TRC / TRA arrays 530 ! 531 #if defined key_diadct 532 ierr = ierr + diadct_alloc () ! 533 #endif 527 534 ! 528 535 IF( lk_mpp ) CALL mpp_sum( ierr )
Note: See TracChangeset
for help on using the changeset viewer.