- Timestamp:
- 2012-11-19T15:13:26+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
r3601 r3602 21 21 !!---------------------------------------------------------------------- 22 22 !!---------------------------------------------------------------------- 23 !! dia_dct : compute the transport through a sec.24 !! dia_dct_init : read namelist.25 !! readsec : read sections description and pathway26 !! removepoints : remove points which are common to 2 procs23 !! dia_dct : Compute the transport through a sec. 24 !! dia_dct_init : Read namelist. 25 !! readsec : Read sections description and pathway 26 !! removepoints : Remove points which are common to 2 procs 27 27 !! transport : Compute transport for each sections 28 !! dia_dct_wri : write tranports results in ascii files29 !! interp : compute Temperature/Salinity/density onU-point or V-point28 !! dia_dct_wri : Write tranports results in ascii files 29 !! interp : Compute temperature/salinity/density at U-point or V-point 30 30 !! 31 31 !!---------------------------------------------------------------------- … … 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 … … 72 73 INTEGER, PARAMETER :: nb_sec_max = 150 73 74 INTEGER, PARAMETER :: nb_point_max = 2000 74 INTEGER, PARAMETER :: nb_type_class = 14 75 INTEGER, PARAMETER :: nb_type_class = 10 76 INTEGER, PARAMETER :: nb_3d_vars = 3 77 INTEGER, PARAMETER :: nb_2d_vars = 2 75 78 INTEGER :: nb_sec 76 79 … … 92 95 INTEGER :: nb_class ! number of boundaries for density classes 93 96 INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section 94 CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! c aracteristics of the class97 CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! characteristics of the class 95 98 REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want) 96 99 zsigp ,&! potential density classes (99 if you don't want) … … 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 … … 190 223 ENDIF 191 224 225 ! Initialise arrays 226 zwork(:) = 0.0 227 zsum(:,:,:) = 0.0 228 192 229 IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN 193 230 WRITE(numout,*) " " … … 208 245 209 246 !Compute transport through section 210 CALL transport(secs(jsec),lldebug )247 CALL transport(secs(jsec),lldebug,jsec) 211 248 212 249 ENDDO … … 214 251 IF( MOD(kt,nn_dctwri)==0 )THEN 215 252 216 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: write at kt = ",kt253 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average transports and write at kt = ",kt 217 254 255 !! divide arrays by nn_dctwri/nn_dct to obtain average 256 transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 257 transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct) 258 259 ! Sum over each class 260 DO jsec=1,nb_sec 261 CALL dia_dct_sum(secs(jsec),jsec) 262 ENDDO 263 218 264 !Sum on all procs 219 265 IF( lk_mpp )THEN … … 233 279 234 280 !nullify transports values after writing 281 transports_3d(:,jsec,:,:)=0. 282 transports_2d(:,jsec,: )=0. 235 283 secs(jsec)%transport(:,:)=0. 236 284 … … 265 313 INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer 266 314 INTEGER :: jsec, jpt ! dummy loop indices 267 ! heat/salt tranport is actived268 315 269 316 INTEGER, DIMENSION(2) :: icoord … … 457 504 !! *** function removepoints 458 505 !! 459 !! ** Purpose :: 460 !! remove points which are common to 2 procs 461 !! 506 !! ** Purpose :: Remove points which are common to 2 procs 462 507 !! 463 508 !---------------------------------------------------------------------------- … … 535 580 END SUBROUTINE removepoints 536 581 537 SUBROUTINE transport(sec,ld_debug )582 SUBROUTINE transport(sec,ld_debug,jsec) 538 583 !!------------------------------------------------------------------------------------------- 539 584 !! *** ROUTINE transport *** 540 585 !! 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. 586 !! Purpose :: Compute the transport for each point in a section 552 587 !! 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 !! 588 !! Method :: Loop over each segment, and each vertical level and add the transport 589 !! Be aware : 590 !! One section is a sum of segments 591 !! One segment is defined by 2 consecutive points in sec%listPoint 592 !! All points of sec%listPoint are positioned on the F-point of the cell 593 !! 594 !! There are two loops: 595 !! loop on the segment between 2 nodes 596 !! loop on the level jk !! 597 !! 598 !! Output :: Arrays containing the volume,density,heat,salt transports for each i 599 !! point in a section, summed over each nn_dct. 561 600 !! 562 601 !!------------------------------------------------------------------------------------------- … … 564 603 TYPE(SECTION),INTENT(INOUT) :: sec 565 604 LOGICAL ,INTENT(IN) :: ld_debug 605 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 566 606 567 607 !! * 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 608 INTEGER :: jk, jseg, jclass, &!loop on level/segment/classes 609 isgnu, isgnv ! 610 REAL(wp) :: zumid, zvmid, &!U/V velocity on a cell segment 611 zumid_ice, zvmid_ice, &!U/V ice velocity 612 zTnorm !transport of velocity through one cell's sides 613 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point 582 614 583 615 TYPE(POINT_SECTION) :: k 584 REAL(wp), POINTER, DIMENSION(:,:):: zsum ! 2D work array585 616 !!-------------------------------------------------------- 586 CALL wrk_alloc( nb_type_class , nb_class_max , zsum )587 617 588 618 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 619 597 620 !---------------------------! … … 670 693 END SELECT 671 694 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) 695 !---------------------------| 696 ! LOOP ON THE LEVEL | 697 !---------------------------| 698 !Sum of the transport on the vertical 699 DO jk=1,mbathy(k%I,k%J) 700 701 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 702 SELECT CASE( sec%direction(jseg) ) 703 CASE(0,1) 704 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 705 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 706 zrhop = interp(k%I,k%J,jk,'V',rhop) 707 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 708 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 709 CASE(2,3) 710 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 711 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 712 zrhop = interp(k%I,k%J,jk,'U',rhop) 713 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 714 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 715 END SELECT 716 717 zfsdep= gdept(k%I,k%J,jk) 718 719 !compute velocity with the correct direction 720 SELECT CASE( sec%direction(jseg) ) 721 CASE(0,1) 722 zumid=0. 723 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 724 CASE(2,3) 725 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 726 zvmid=0. 727 END SELECT 728 729 !zTnorm=transport through one cell; 730 !velocity* cell's length * cell's thickness 731 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 732 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 740 733 741 734 #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 ENDIF735 !add transport due to free surface 736 IF( jk==1 )THEN 737 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 738 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 739 ENDIF 747 740 #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 741 !COMPUTE TRANSPORT 742 743 transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 744 745 IF ( sec%llstrpond ) THEN 746 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp 747 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001 792 748 ENDIF 793 749 794 ENDDO !end of loop on the density classes750 ENDDO !end of loop on the level 795 751 796 752 #if defined key_lim2 || defined key_lim3 … … 816 772 zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 817 773 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 774 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* & 775 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 776 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 777 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 778 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* & 779 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 842 780 843 781 ENDIF !end of ice case … … 846 784 ENDDO !end of loop on the segment 847 785 848 849 ELSE !if sec%nb_point =0 850 zsum(1:2,:)=0. 851 IF (sec%llstrpond) zsum(3:10,:)=0. 852 zsum( 11:14,:)=0. 853 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-6 860 sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6 861 IF( sec%llstrpond ) THEN 862 IF( zsum(1,jclass) .NE. 0._wp ) THEN 863 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 ENDIF 868 IF( zsum(2,jclass) .NE. 0._wp )THEN 869 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 ENDIF 874 ELSE 875 sec%transport( 3,jclass) = 0._wp 876 sec%transport( 4,jclass) = 0._wp 877 sec%transport( 5,jclass) = 0._wp 878 sec%transport( 6,jclass) = 0._wp 879 sec%transport( 7,jclass) = 0._wp 880 sec%transport( 8,jclass) = 0._wp 881 sec%transport(10,jclass) = 0._wp 882 ENDIF 883 ENDDO 884 885 IF( sec%ll_ice_section ) THEN 886 sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6 887 sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6 888 sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6 889 sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6 890 ENDIF 891 892 CALL wrk_dealloc( nb_type_class , nb_class_max , zsum ) 786 ENDIF !end of sec%nb_point =0 case 893 787 ! 894 788 END SUBROUTINE transport 789 790 SUBROUTINE dia_dct_sum(sec,jsec) 791 !!------------------------------------------------------------- 792 !! Purpose: Average the transport over nn_dctwri time steps 793 !! and sum over the density/salinity/temperature/depth classes 794 !! 795 !! Method: Sum over relevant grid cells to obtain values 796 !! for each class 797 !! There are several loops: 798 !! loop on the segment between 2 nodes 799 !! loop on the level jk 800 !! loop on the density/temperature/salinity/level classes 801 !! test on the density/temperature/salinity/level 802 !! 803 !! Note: Transport through a given section is equal to the sum of transports 804 !! computed on each proc. 805 !! On each proc,transport is equal to the sum of transport computed through 806 !! segments linking each point of sec%listPoint with the next one. 807 !! 808 !!------------------------------------------------------------- 809 !! * arguments 810 TYPE(SECTION),INTENT(INOUT) :: sec 811 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 812 813 TYPE(POINT_SECTION) :: k 814 INTEGER :: jk,jseg,jclass ! dummy variables for looping on level/segment/classes 815 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 816 !!------------------------------------------------------------- 817 818 !! Sum the relevant segments to obtain values for each class 819 IF(sec%nb_point .NE. 0)THEN 820 821 !--------------------------------------! 822 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 823 !--------------------------------------! 824 DO jseg=1,MAX(sec%nb_point-1,0) 825 826 !------------------------------------------------------------------------------------------- 827 ! Select the appropriate coordinate for computing the velocity of the segment 828 ! 829 ! CASE(0) Case (2) 830 ! ------- -------- 831 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 832 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 833 ! | 834 ! | 835 ! | 836 ! Case (3) U(i,j) 837 ! -------- | 838 ! | 839 ! listPoint(jseg+1) F(i,j+1) | 840 ! | | 841 ! | | 842 ! | listPoint(jseg+1) F(i,j-1) 843 ! | 844 ! | 845 ! U(i,j+1) 846 ! | Case(1) 847 ! | ------ 848 ! | 849 ! | listPoint(jseg+1) listPoint(jseg) 850 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 851 ! listPoint(jseg) F(i,j) 852 ! 853 !------------------------------------------------------------------------------------------- 854 855 SELECT CASE( sec%direction(jseg) ) 856 CASE(0) ; k = sec%listPoint(jseg) 857 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 858 CASE(2) ; k = sec%listPoint(jseg) 859 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 860 END SELECT 861 862 !---------------------------| 863 ! LOOP ON THE LEVEL | 864 !---------------------------| 865 !Sum of the transport on the vertical 866 DO jk=1,mbathy(k%I,k%J) 867 868 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 869 SELECT CASE( sec%direction(jseg) ) 870 CASE(0,1) 871 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 872 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 873 zrhop = interp(k%I,k%J,jk,'V',rhop) 874 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 875 876 CASE(2,3) 877 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 878 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 879 zrhop = interp(k%I,k%J,jk,'U',rhop) 880 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 881 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 882 END SELECT 883 884 zfsdep= gdept(k%I,k%J,jk) 885 886 !------------------------------- 887 ! LOOP ON THE DENSITY CLASSES | 888 !------------------------------- 889 !The computation is made for each density/temperature/salinity/depth class 890 DO jclass=1,MAX(1,sec%nb_class-1) 891 892 !----------------------------------------------! 893 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 894 !----------------------------------------------! 895 896 IF ( ( & 897 ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 898 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 899 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 900 901 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 902 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 903 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 904 905 ((( zsn .GT. sec%zsal(jclass)) .AND. & 906 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 907 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 908 909 ((( ztn .GE. sec%ztem(jclass)) .AND. & 910 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 911 ( sec%ztem(jclass) .EQ.99.)) .AND. & 912 913 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 914 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 915 ( sec%zlay(jclass) .EQ. 99. )) & 916 )) THEN 917 918 !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 919 !---------------------------------------------------------------------------- 920 IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 921 sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 922 ELSE 923 sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 924 ENDIF 925 IF( sec%llstrpond )THEN 926 927 IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN 928 sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 929 ELSE 930 sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 931 ENDIF 932 933 IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN 934 sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 935 ELSE 936 sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 937 ENDIF 938 939 ELSE 940 sec%transport( 3,jclass) = 0._wp 941 sec%transport( 4,jclass) = 0._wp 942 sec%transport( 5,jclass) = 0._wp 943 sec%transport( 6,jclass) = 0._wp 944 ENDIF 945 946 ENDIF ! end of test if point is in class 947 948 ENDDO ! end of loop on the classes 949 950 ENDDO ! loop over jk 951 952 #if defined key_lim2 || defined key_lim3 953 954 !ICE CASE 955 IF( sec%ll_ice_section )THEN 956 957 IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 958 sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6 959 ELSE 960 sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6 961 ENDIF 962 963 IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 964 sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6 965 ELSE 966 sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6 967 ENDIF 968 969 ENDIF !end of ice case 970 #endif 971 972 ENDDO !end of loop on the segment 973 974 ELSE !if sec%nb_point =0 975 sec%transport(1:2,:)=0. 976 IF (sec%llstrpond) sec%transport(3:6,:)=0. 977 IF (sec%ll_ice_section) sec%transport(7:10,:)=0. 978 ENDIF !end of sec%nb_point =0 case 979 980 END SUBROUTINE dia_dct_sum 895 981 896 982 SUBROUTINE dia_dct_wri(kt,ksec,sec) … … 905 991 !! 906 992 !! 2. Write heat transports in "heat_transport" 907 !! Unit: Peta W : area * Velocity * T * rh au * Cp / 1.e15993 !! Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15 908 994 !! 909 995 !! 3. Write salt transports in "salt_transport" 910 !! Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6996 !! Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9 911 997 !! 912 998 !!------------------------------------------------------------- … … 917 1003 918 1004 !!local declarations 919 INTEGER :: jcl ,ji! Dummy loop1005 INTEGER :: jclass ! Dummy loop 920 1006 CHARACTER(len=2) :: classe ! Classname 921 1007 REAL(wp) :: zbnd1,zbnd2 ! Class bounds 922 1008 REAL(wp) :: zslope ! section's slope coeff 923 1009 ! 924 REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace1010 REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 925 1011 !!------------------------------------------------------------- 926 CALL wrk_alloc(nb_type_class , zsumclass )927 928 zsumclass (:)=0._wp1012 CALL wrk_alloc(nb_type_class , zsumclasses ) 1013 1014 zsumclasses(:)=0._wp 929 1015 zslope = sec%slopeSection 930 1016 931 1017 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) 1018 DO jclass=1,MAX(1,sec%nb_class-1) 1019 936 1020 classe = 'N ' 937 1021 zbnd1 = 0._wp 938 1022 zbnd2 = 0._wp 939 zsumclass (1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)1023 zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass) 940 1024 941 1025 942 1026 !insitu density classes transports 943 IF( ( sec%zsigi(jcl ) .NE. 99._wp ) .AND. &944 ( sec%zsigi(jcl +1) .NE. 99._wp ) )THEN1027 IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. & 1028 ( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN 945 1029 classe = 'DI ' 946 zbnd1 = sec%zsigi(jcl )947 zbnd2 = sec%zsigi(jcl +1)1030 zbnd1 = sec%zsigi(jclass) 1031 zbnd2 = sec%zsigi(jclass+1) 948 1032 ENDIF 949 1033 !potential density classes transports 950 IF( ( sec%zsigp(jcl ) .NE. 99._wp ) .AND. &951 ( sec%zsigp(jcl +1) .NE. 99._wp ) )THEN1034 IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. & 1035 ( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN 952 1036 classe = 'DP ' 953 zbnd1 = sec%zsigp(jcl )954 zbnd2 = sec%zsigp(jcl +1)1037 zbnd1 = sec%zsigp(jclass) 1038 zbnd2 = sec%zsigp(jclass+1) 955 1039 ENDIF 956 1040 !depth classes transports 957 IF( ( sec%zlay(jcl ) .NE. 99._wp ) .AND. &958 ( sec%zlay(jcl +1) .NE. 99._wp ) )THEN1041 IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. & 1042 ( sec%zlay(jclass+1) .NE. 99._wp ) )THEN 959 1043 classe = 'Z ' 960 zbnd1 = sec%zlay(jcl )961 zbnd2 = sec%zlay(jcl +1)1044 zbnd1 = sec%zlay(jclass) 1045 zbnd2 = sec%zlay(jclass+1) 962 1046 ENDIF 963 1047 !salinity classes transports 964 IF( ( sec%zsal(jcl ) .NE. 99._wp ) .AND. &965 ( sec%zsal(jcl +1) .NE. 99._wp ) )THEN1048 IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. & 1049 ( sec%zsal(jclass+1) .NE. 99._wp ) )THEN 966 1050 classe = 'S ' 967 zbnd1 = sec%zsal(jcl )968 zbnd2 = sec%zsal(jcl +1)1051 zbnd1 = sec%zsal(jclass) 1052 zbnd2 = sec%zsal(jclass+1) 969 1053 ENDIF 970 1054 !temperature classes transports 971 IF( ( sec%ztem(jcl ) .NE. 99._wp ) .AND. &972 ( sec%ztem(jcl +1) .NE. 99._wp ) ) THEN1055 IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. & 1056 ( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN 973 1057 classe = 'T ' 974 zbnd1 = sec%ztem(jcl )975 zbnd2 = sec%ztem(jcl +1)1058 zbnd1 = sec%ztem(jclass) 1059 zbnd2 = sec%ztem(jclass+1) 976 1060 ENDIF 977 1061 978 1062 !write volume transport per class 979 1063 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)1064 jclass,classe,zbnd1,zbnd2,& 1065 sec%transport(1,jclass),sec%transport(2,jclass), & 1066 sec%transport(1,jclass)+sec%transport(2,jclass) 983 1067 984 1068 IF( sec%llstrpond )THEN … … 986 1070 !write heat transport per class: 987 1071 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.e151072 jclass,classe,zbnd1,zbnd2,& 1073 sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, & 1074 ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15 991 1075 !write salt transport per class 992 1076 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.e91077 jclass,classe,zbnd1,zbnd2,& 1078 sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,& 1079 (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9 996 1080 ENDIF 997 1081 … … 1000 1084 zbnd1 = 0._wp 1001 1085 zbnd2 = 0._wp 1002 jcl =01086 jclass=0 1003 1087 1004 1088 !write total volume transport 1005 1089 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 1006 jcl ,"total",zbnd1,zbnd2,&1007 zsumclass (1),zsumclass(2),zsumclass(1)+zsumclass(2)1090 jclass,"total",zbnd1,zbnd2,& 1091 zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2) 1008 1092 1009 1093 IF( sec%llstrpond )THEN … … 1011 1095 !write total heat transport 1012 1096 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.e151097 jclass,"total",zbnd1,zbnd2,& 1098 zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,& 1099 (zsumclasses(3)+zsumclasses(4) )*1.e-15 1016 1100 !write total salt transport 1017 1101 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.e91102 jclass,"total",zbnd1,zbnd2,& 1103 zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,& 1104 (zsumclasses(5)+zsumclasses(6))*1.e-9 1021 1105 ENDIF 1022 1106 … … 1025 1109 !write total ice volume transport 1026 1110 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1027 jcl ,"ice_vol",zbnd1,zbnd2,&1028 sec%transport( 9,1),sec%transport(10,1),&1029 sec%transport( 9,1)+sec%transport(10,1)1111 jclass,"ice_vol",zbnd1,zbnd2,& 1112 sec%transport(7,1),sec%transport(8,1),& 1113 sec%transport(7,1)+sec%transport(8,1) 1030 1114 !write total ice surface transport 1031 1115 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1032 jcl ,"ice_surf",zbnd1,zbnd2,&1033 sec%transport( 11,1),sec%transport(12,1), &1034 sec%transport( 11,1)+sec%transport(12,1)1116 jclass,"ice_surf",zbnd1,zbnd2,& 1117 sec%transport(9,1),sec%transport(10,1), & 1118 sec%transport(9,1)+sec%transport(10,1) 1035 1119 ENDIF 1036 1120 … … 1038 1122 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 1039 1123 1040 CALL wrk_dealloc(nb_type_class , zsumclass )1124 CALL wrk_dealloc(nb_type_class , zsumclasses ) 1041 1125 END SUBROUTINE dia_dct_wri 1042 1126 … … 1044 1128 !!---------------------------------------------------------------------- 1045 1129 !! 1046 !! Purpose: compute Temperature/Salinity/density at U-point or V-point1130 !! Purpose: compute temperature/salinity/density at U-point or V-point 1047 1131 !! -------- 1048 1132 !! … … 1053 1137 !! 1054 1138 !! 1055 !! | I | I+1 | Z= Temperature/Salinity/density at U-poinT1139 !! | I | I+1 | Z=temperature/salinity/density at U-poinT 1056 1140 !! | | | 1057 !! ---------------------------------------- 1. Veritcal einterpolation: compute zbis1141 !! ---------------------------------------- 1. Veritcal interpolation: compute zbis 1058 1142 !! | | | interpolation between ptab(I,J,K) and ptab(I,J,K+1) 1059 1143 !! | | | zbis = … … 1136 1220 zdep2 = fsdept(ii2,ij2,kk) - zdepu 1137 1221 1138 ! weights1222 ! weights 1139 1223 zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) ) 1140 1224 zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) ) … … 1163 1247 1164 1248 IF( ze3t >= 0. )THEN 1165 ! zbis1249 ! zbis 1166 1250 zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 1167 1251 ! result 1168 1252 interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 1169 1253 ELSE 1170 ! zbis1254 ! zbis 1171 1255 zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 1172 1256 ! result … … 1195 1279 END SUBROUTINE dia_dct_init 1196 1280 1197 SUBROUTINE dia_dct( kt ) 1198 INTEGER, INTENT( in ) :: kt! ocean time-step index1281 SUBROUTINE dia_dct( kt ) ! Dummy routine 1282 INTEGER, INTENT( in ) :: kt ! ocean time-step index 1199 1283 WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt 1200 1284 END SUBROUTINE dia_dct
Note: See TracChangeset
for help on using the changeset viewer.
