MODULE diadct !!===================================================================== !! *** MODULE diadct *** !! Ocean diagnostics: Compute the transport trough a sec. !!=============================================================== !! History : !! !! !! !!---------------------------------------------------------------------- #if defined key_diadct !!---------------------------------------------------------------------- !! 'key_diadct' : !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dia_dct : Compute the transport trough a sec. !! dia_dct_init : !! readsec : !! removepoints : !! transport : !! dia_dct_wri : !! interp : !! !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE in_out_manager ! I/O manager USE daymod ! calendar USE dianam ! build name of file USE lib_mpp ! distributed memory computing library #if defined key_ice_lim2 || defined key_ice_lim3 USE ice #endif IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC dia_dct ! routine called by step.F90 PUBLIC dia_dct_init! routine called by opa.F90 PRIVATE readsec PRIVATE removepoints PRIVATE transport PRIVATE dia_dct_wri #include "domzgr_substitute.h90" !! * Shared module variables LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag !! * Module variables INTEGER :: nn_dct = 1 ! Frequency of computation INTEGER :: nn_dctwri = 1 ! Frequency of output INTEGER :: nn_secdebug = 0 ! Number of the section to debug INTEGER, PARAMETER :: nb_class_max = 10 INTEGER, PARAMETER :: nb_sec_max = 150 INTEGER, PARAMETER :: nb_point_max = 2000 INTEGER, PARAMETER :: nb_type_class = 14 INTEGER :: nb_sec TYPE POINT_SECTION INTEGER :: I,J END TYPE POINT_SECTION TYPE COORD_SECTION REAL(wp) :: lon,lat END TYPE COORD_SECTION TYPE SECTION CHARACTER(len=60) :: name ! name of the sec LOGICAL :: llstrpond ! true if you want the computation of salinity and ! temperature balanced by the transport LOGICAL :: ll_ice_section ! icesurf and icevol computation LOGICAL :: ll_date_line ! = T if the section crosses the date-line TYPE(COORD_SECTION), DIMENSION(2) :: coordSec ! longitude and latitude of the extremities of the sec INTEGER :: nb_class ! number of boundaries for density classes INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! caracteristics of the class REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want) zsigp ,&! potential density classes (99 if you don't want) zsal ,&! salinity classes (99 if you don't want) ztem ,&! temperature classes(99 if you don't want) zlay ! level classes (99 if you don't want) REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output REAL(wp) :: slopeSection ! coeff directeur de la section INTEGER :: nb_point ! nb de points de la section TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of point in the section END TYPE SECTION TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections CONTAINS SUBROUTINE dia_dct_init !!--------------------------------------------------------------------- !! *** ROUTINE diadct *** !! !! ** Purpose: Read the namelist parametres and the !! file containing the sections definitions !! !! ** Input : Read a binary file generated by the tools "????" !! !! ** Output: secs array !! !! History : !! ! 03-09 (JM. Lellouche and E. Durand) mpprsum modified !! ! 04-01 (C.Derval) bugfix !! ! 04-2007 (C.Bricaud) bugfix !! ! mise à zéro de secs(jsec)%transport(:,:) apres chaque ecriture !!--------------------------------------------------------------------- NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug !read namelist REWIND( numnam ) READ ( numnam, namdct ) IF( lwp ) THEN WRITE(numout,*) " " WRITE(numout,*) "diadct_init: compute transports through sections " WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" WRITE(numout,*) " Frequency of computation: nn_dct = ",nn_dct WRITE(numout,*) " Frequency of write: nn_dctwri = ",nn_dctwri IF ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN WRITE(numout,*)" Debug section number: ", nn_secdebug ELSE IF ( nn_secdebug == 0 )THEN ; WRITE(numout,*)" No section to debug" ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)" Debug all sections" ELSE ; WRITE(numout,*)" Wrong value for nn_secdebug : ",nn_secdebug ENDIF IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0) & & CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' ) !open output file CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) CALL ctl_opn( numdct_heat, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) CALL ctl_opn( numdct_salt, 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) ENDIF !Read section_ijglobal.diadct CALL readsec END SUBROUTINE dia_dct_init SUBROUTINE dia_dct(kt) !!--------------------------------------------------------------------- !! *** ROUTINE diadct *** !! !! ** Purpose: Compute sections tranport and write it in numdct file !!--------------------------------------------------------------------- !! * Arguments INTEGER,INTENT(IN) ::kt !! * Local variables INTEGER :: jsec, &!loop on sections iost !error for opening fileout LOGICAL :: lldebug =.FALSE. !debug a section CHARACTER(len=160) :: clfileout !fileout name !!--------------------------------------------------------------------- IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN WRITE(numout,*) " " WRITE(numout,*) "diadct: compute transport" WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~" WRITE(numout,*) "nb_sec = ",nb_sec ENDIF ! Compute transport and write only at nn_dctwri IF( MOD(kt,nn_dct)==0 ) THEN DO jsec=1,nb_sec !debug this section computing ? lldebug=.FALSE. ! IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 ) lldebug=.TRUE. !Compute transport through section CALL transport(secs(jsec),lldebug) IF( MOD(kt,nn_dctwri)==0 )THEN IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: write at kt = ",kt !Write the transport IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec)) !nullify transports values after writing secs(jsec)%transport(:,:)=0. ENDIF ENDDO ENDIF END SUBROUTINE dia_dct SUBROUTINE readsec !!--------------------------------------------------------------------- !! *** ROUTINE readsec *** !! !! ** Purpose: Read the sections definitions !!--------------------------------------------------------------------- !! * Local variables INTEGER :: iptglo , iptloc ! Global and local number of points for a section INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer INTEGER :: jsec, jpt ! dummy loop indices INTEGER, DIMENSION(2) :: icoord CHARACTER(len=160) :: clname !filename CHARACTER(len=200) :: cltmp TYPE(POINT_SECTION),DIMENSION(nb_point_max) ::coordtemp !contains listpoints coordinates !read in the file INTEGER,DIMENSION(nb_point_max) ::directemp !contains listpoints directions !read in the files LOGICAL :: llbon ,&!local logical lldebug !debug the section !!------------------------------------------------------------------------------------- !open input file !--------------- CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) !--------------- !Read input file !--------------- DO jsec=1,nb_sec_max !loop on the nb_sec sections IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) & & WRITE(numout,*)'debuging for section number: ',jsec !initialization !--------------- secs(jsec)%name='' secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE. secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0 secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp secs(jsec)%zlay = 99._wp secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0 !read section's number / name / computing choices / classes / slopeSection / points number !----------------------------------------------------------------------------------------- READ(numdct_in,iostat=iost)isec IF (iost .NE. 0 )EXIT !end of file WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' et jsec= ',jsec IF( jsec .NE. isec ) CALL ctl_stop( cltmp ) IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec READ(numdct_in)secs(jsec)%name READ(numdct_in)secs(jsec)%llstrpond READ(numdct_in)secs(jsec)%ll_ice_section READ(numdct_in)secs(jsec)%ll_date_line READ(numdct_in)secs(jsec)%coordSec READ(numdct_in)secs(jsec)%nb_class READ(numdct_in)secs(jsec)%zsigi READ(numdct_in)secs(jsec)%zsigp READ(numdct_in)secs(jsec)%zsal READ(numdct_in)secs(jsec)%ztem READ(numdct_in)secs(jsec)%zlay READ(numdct_in)secs(jsec)%slopeSection READ(numdct_in)iptglo !debug !----- IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name) WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line WRITE(numout,*) " Slope section : ",secs(jsec)%slopeSection WRITE(numout,*) " Number of points in the section: ",iptglo WRITE(numout,*) " Number of classes ",secs(jsec)%nb_class WRITE(numout,'(A40,10(f8.3,1X))')" Insitu density classes : ",secs(jsec)%zsigi WRITE(numout,'(A40,10(f8.3,1X))')" Potential density classes : ",secs(jsec)%zsigp WRITE(numout,'(A40,10(f8.3,1X))')" Salinity classes : ",secs(jsec)%zsal WRITE(numout,'(A40,10(f8.3,1X))')" Temperature classes : ",secs(jsec)%ztem WRITE(numout,'(A40,10(f8.3,1X))')" Depth classes : ",secs(jsec)%zlay ENDIF IF( iptglo .NE. 0 )THEN !read points'coordinates and directions !-------------------------------------- coordtemp(:) = POINT_SECTION(0,0) !list of points read directemp(:) = 0 !value of directions of each points DO jpt=1,iptglo READ(numdct_in)i1,i2 coordtemp(jpt)%I = i1 coordtemp(jpt)%J = i2 ENDDO READ(numdct_in)directemp(1:iptglo) !debug !----- IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN WRITE(numout,*)" List of points in global domain:" DO jpt=1,iptglo WRITE(numout,*)' # I J ',jpt,coordtemp(jpt) ENDDO ENDIF !?IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN ! DO jpt=1,iptglo ! WRITE(numout,*)' # I J ',jpt,coordtemp(jpt) ! ENDDO !ENDIF !Now each proc selects only points that are in its domain: !-------------------------------------------------------- iptloc = 0 !initialize number of points selected DO jpt=1,iptglo !loop on listpoint read in the file iiglo=coordtemp(jpt)%I ! global coordinates of the point ijglo=coordtemp(jpt)%J ! " IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2 iiloc=iiglo-jpizoom+1-nimpp+1 ! local coordinates of the point ijloc=ijglo-jpjzoom+1-njmpp+1 ! " !verify if the point is on the local domain:(1,nlei)*(1,nlej) IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. & ijloc .GE. 1 .AND. ijloc .LE. nlej )THEN iptloc = iptloc + 1 ! count local points secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates secs(jsec)%direction(iptloc) = directemp(jpt) ! store local direction ENDIF ENDDO secs(jsec)%nb_point=iptloc !store number of section's points !debug !----- IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN WRITE(numout,*)" List of points selected by the proc:" DO jpt = 1,iptloc iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 WRITE(numout,*)' # I J : ',iiglo,ijglo ENDDO ENDIF !remove redundant points between processors !look NEMO documentation for more explanations ! lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE. lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) ) lldebug = .TRUE. IF( iptloc .NE. 0 )THEN CALL removepoints(secs(jsec),'I','top_list',lldebug) CALL removepoints(secs(jsec),'I','bot_list',lldebug) CALL removepoints(secs(jsec),'J','top_list',lldebug) CALL removepoints(secs(jsec),'J','bot_list',lldebug) ENDIF !debug !----- IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN WRITE(numout,*)" List of points after removepoints:" DO jpt = 1,iptloc iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 WRITE(numout,*)' # I J : ',iiglo,ijglo ENDDO ENDIF ELSE ! iptglo = 0 IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )& WRITE(numout,*)' Not points for this section.' ENDIF ENDDO !end of the loop on jsec nb_sec = jsec-1 !number of section read in the file END SUBROUTINE readsec SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug) !!--------------------------------------------------------------------------- !! *** function removepoints !! !---------------------------------------------------------------------------- !! * arguments TYPE(SECTION),INTENT(INOUT) :: sec CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J' CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list' LOGICAL,INTENT(IN) :: ld_debug !! * Local variables INTEGER :: iextr ,& !extremity of listpoint that we verify iind ,& !coord of listpoint that we verify itest ,& !indice value of the side of the domain !where points could be redundant isgn ,& !way of course in listpoint istart,iend !first and last points selected in listpoint INTEGER :: jpoint =0 !loop on list points INTEGER,DIMENSION(nb_point_max) :: idirec !contains temporare sec%direction INTEGER,DIMENSION(2,nb_point_max) :: icoord !contains temporare sec%listpoint !---------------------------------------------------------------------------- IF( ld_debug )WRITE(numout,*)' -------------------------' IF( ld_debug )WRITE(numout,*)' removepoints in listpoint' !iextr=extremity of list_point that we verify IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1 ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1 ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr") ENDIF !which coordinate shall we verify ? IF ( cdind=='I' )THEN ; itest=nlei ; iind=1 ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2 ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind") ENDIF IF( ld_debug )THEN WRITE(numout,*)' case: coord/list extr/domain side' WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest WRITE(numout,*)' Actual number of points: ',sec%nb_point ENDIF icoord(1,1:nb_point_max) = sec%listPoint%I icoord(2,1:nb_point_max) = sec%listPoint%J idirec = sec%direction sec%listPoint = POINT_SECTION(0,0) sec%direction = 0 jpoint=iextr+isgn DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. & icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest) jpoint=jpoint+isgn ENDDO IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point ELSE ; istart=1 ; iend=jpoint+1 ENDIF sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend) sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend) sec%direction(1:1+iend-istart) = idirec(istart:iend) sec%nb_point = iend-istart+1 IF( ld_debug )THEN WRITE(numout,*)' Number of points after removepoints :',sec%nb_point WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point) ENDIF END SUBROUTINE removepoints SUBROUTINE transport(sec,ld_debug) !!--------------------------------------------------------------------- !! *** ROUTINE transport *** !! !! ** Purpose : Compute the transport trough a sec !! !! ** Method :Transport through a given section is equal to the sum of transports !! computed on each proc. !! On each proc,transport is equal to the sum of transport computed through !! segments linking each points of sec%listPoint with the next one. !! There are several loops: !! loop on the density/temperature/salinity/level classes !! loop on the segment between 2 nodes !! loop on the level jk !! test on the density/temperature/salinity/level !! !! ** Output: sec%transport: transport in the 2 direction and temperature, !! salinity, density meaned by the transport !! ** History : !! original : 02/99 (Y Drillet) !! addition : 10/01 (Y Drillet, R Bourdalle Badie) !! : 10/05 (M Laborie) F90 !! addition : 04/07 (G Garric) Ice sections !! bugfix : 04/07 (C Bricaud) test on sec%nb_point !! initialisation of ztransp1,ztransp2,... !!------------------------------------------------------------------- !! * Arguments TYPE(SECTION),INTENT(INOUT) :: sec LOGICAL ,INTENT(IN) :: ld_debug !! * Local variables INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes iptegrid, &!=0 if section is constant along i/j isgnu , isgnv ! INTEGER :: ii, ij ! local integer REAL(wp):: zumid , zvmid ,&!U/V velocity on a cell segment zumid_ice , zvmid_ice ,&!U/V ice velocity zTnorm ,&!transport of velocity through one cell's sides ztransp1 , ztransp2 ,&!total transport in sens 1 and 2 ztemp1 , ztemp2 ,&!temperature transport " zrhoi1 , zrhoi2 ,&!mass transport " zrhop1 , zrhop2 ,&!mass transport " zsal1 , zsal2 ,&!salinity transport " zice_vol_pos , zice_vol_neg ,&!volumic ice transport " zice_surf_pos, zice_surf_neg !surfacic ice transport " REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point REAL(wp) :: aj,i0,j0,i1,j1,i,j TYPE(POINT_SECTION) :: k INTEGER,DIMENSION(1):: ish INTEGER,DIMENSION(2):: ish2 REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork ! temporary vector for mpp_sum !!-------------------------------------------------------- IF( ld_debug )WRITE(numout,*)' Compute transport' !----------------! ! INITIALIZATION ! !----------------! zsum = 0._wp zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp zice_vol_pos = 0._wp ; zice_vol_neg = 0._wp !compute iptegrid: iptegrid=0 if section is constant along i or j iptegrid = 0 jseg = 1 DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 ) IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1 jseg = jseg+1 ENDDO IF(lk_mpp) CALL mpp_sum(iptegrid) !---------------------------! ! COMPUTE TRANSPORT ! !---------------------------! IF(sec%nb_point .NE. 0)THEN !---------------------------------------------------------------------------------------------------- !Compute sign for velocities: ! !convention: ! non horizontal section: sens + is toward left hand of section ! horizontal section: sens + is toward north of section ! ! ! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0 ! ---------------- ----------------- --------------- -------------- ! ! isgnv=1 sens + ! ______ _____ ______(i1,j1) (i1,j1) !(i0,j0) | //| | | sens + ! | isgnu=1 // | |isgnu=1 |isgnu=1 /|\ ! |_______ // ______| \\ | ---\ | ! | | isgnv=-1 \\ | | ---/ sens + ____________(i1,j1) ! | | __\\| | (i0,j0) ! |(i1,j1) |(i0,j0) sens + | isgnv=1 ! (i0,j0) !---------------------------------------------------------------------------------------------------- isgnu = 1 IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1 ELSE ; isgnv = 1 ENDIF IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv !--------------------------------------! ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! !--------------------------------------! DO jseg=1,MAX(sec%nb_point-1,0) !------------------------------------------ ! Select good coordinate to have velocity of the segment ! ! |U CASE(3) ! CASE(0) | j+1 ! ____V_____|F____V_CASE(1)__ ! | ! | j ! T |U CASE(2) ! i | i+1 !------------------------------------------- SELECT CASE( sec%direction(jseg) ) CASE(0) ; k = sec%listPoint(jseg) CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) CASE(2) ; k = sec%listPoint(jseg) CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) END SELECT !------------------------------- ! LOOP ON THE DENSITY CLASSES | !------------------------------- !The computation is made for each density class DO jclass=1,MAX(1,sec%nb_class-1) ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp !---------------------------| ! LOOP ON THE LEVEL | !---------------------------| !Sum of the transport on the vertical DO jk=1,jpk ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point SELECT CASE( sec%direction(jseg) ) CASE(0,1) ztn = interp(k%I,k%J,jk,'V',tn) zsn = interp(k%I,k%J,jk,'V',sn) zrhop = interp(k%I,k%J,jk,'V',rhop) zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) CASE(2,3) ztn = interp(k%I,k%J,jk,'U',tn) zsn = interp(k%I,k%J,jk,'U',sn) zrhop = interp(k%I,k%J,jk,'U',rhop) zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) END SELECT zfsdep= gdept(k%I,k%J,jk) !----------------------------------------------! !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! !----------------------------------------------! IF ( ( ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & ( sec%zsigp(jclass) .EQ. 99.)) .AND. & ((( zrhoi .GE. (sec%zsigi(jclass) )) .AND. & ( zrhoi .LE. (sec%zsigi(jclass+1)))) .OR. & ( sec%zsigi(jclass) .EQ. 99.)) .AND. & ((( zsn .GT. sec%zsal(jclass)) .AND. & ( zsn .LE. sec%zsal(jclass+1))) .OR. & ( sec%zsal(jclass) .EQ. 99.)) .AND. & ((( ztn .GE. sec%ztem(jclass)) .AND. & ( ztn .LE. sec%ztem(jclass+1))) .OR. & ( sec%ztem(jclass) .EQ.99.)) .AND. & ((( zfsdep .GE. sec%zlay(jclass)) .AND. & ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & ( sec%zlay(jclass) .EQ. 99. )))) THEN !compute velocity with the good sens SELECT CASE( sec%direction(jseg) ) CASE(0,1) zumid=0. zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) CASE(2,3) zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) zvmid=0. END SELECT !velocity* cell's lenght * cell's thickness zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) #if ! defined key_vvl !add transport due to free surface IF( jk==1 )THEN zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) ENDIF #endif !COMPUTE TRANSPORT !zTnorm=transport through one cell for one class !ztransp1 or ztransp2=transport through one cell i ! for one class for one direction IF( zTnorm .GE. 0 )THEN ztransp1=zTnorm+ztransp1 IF ( sec%llstrpond ) THEN ztemp1 = ztemp1 + zTnorm * ztn zsal1 = zsal1 + zTnorm * zsn zrhoi1 = zrhoi1 + zTnorm * zrhoi zrhop1 = zrhop1 + zTnorm * zrhop ENDIF ELSE ztransp2=(zTnorm)+ztransp2 IF ( sec%llstrpond ) THEN ztemp2 = ztemp2 + zTnorm * ztn zsal2 = zsal2 + zTnorm * zsn zrhoi2 = zrhoi2 + zTnorm * zrhoi zrhop2 = zrhop2 + zTnorm * zrhop ENDIF ENDIF ENDIF ! end of density test ENDDO!end of loop on the level !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE DIRECTIONS !--------------------------------------------------- zsum(1,jclass) = zsum(1,jclass)+ztransp1 zsum(2,jclass) = zsum(2,jclass)+ztransp2 IF( sec%llstrpond )THEN zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1 zsum(10,jclass) = zsum(10,jclass)+zsal2 ENDIF ENDDO !end of loop on the density classes #if defined key_ice_lim !ICE CASE !------------ IF( sec%ll_ice_section )THEN SELECT CASE (sec%direction(jseg)) CASE(0) zumid_ice = 0 zvmid_ice = isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) CASE(1) zumid_ice = 0 zvmid_ice = isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) CASE(2) zvmid_ice = 0 zumid_ice = isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) CASE(3) zvmid_ice = 0 zumid_ice = isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) END SELECT zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) IF( zTnorm .GE. 0)THEN zice_vol_pos = (zTnorm)* & (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & +zice_vol_pos zice_surf_pos = (zTnorm)* & (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & +zice_surf_pos ELSE zice_vol_neg=(zTnorm)* & (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & +zice_vol_neg zice_surf_neg=(zTnorm)* & (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & +zice_surf_neg ENDIF zsum(11,1) = zsum(11,1)+zice_vol_pos zsum(12,1) = zsum(12,1)+zice_vol_neg zsum(13,1) = zsum(13,1)+zice_surf_pos zsum(14,1) = zsum(14,1)+zice_surf_neg ENDIF !end of ice case #endif ENDDO !end of loop on the segment ELSE !if nb_inmesh=0 zsum(1:2,:)=0. IF (sec%llstrpond) zsum(3:10,:)=0. zsum( 11:14,:)=0. ENDIF !end of nb_inmesh=0 case !------------------! !SUM ON ALL PROCS ! !------------------! IF( lk_mpp )THEN ish(1) = nb_type_class*nb_class_max ; ish2(1)=nb_type_class ; ish2(2)=nb_class_max zwork(:)= RESHAPE(zsum(:,:), ish ) CALL mpp_sum(zwork, ish(1)) zsum(:,:)= RESHAPE(zwork,ish2) ENDIF !-------------------------------| !FINISH COMPUTING TRANSPORTS | !-------------------------------| DO jclass=1,MAX(1,sec%nb_class-1) sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6 sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6 IF( sec%llstrpond ) THEN IF( zsum(1,jclass) .NE. 0 ) THEN sec%transport(3,jclass)=sec%transport(3,jclass)+zsum(3,jclass)/zsum(1,jclass) sec%transport(5,jclass)=sec%transport(5,jclass)+zsum(5,jclass)/zsum(1,jclass) sec%transport(7,jclass)=sec%transport(7,jclass)+zsum(7,jclass)/zsum(1,jclass) sec%transport(9,jclass)=sec%transport(9,jclass)+zsum(9,jclass)/zsum(1,jclass) ELSE sec%transport(3,jclass)=0. sec%transport(5,jclass)=0. sec%transport(7,jclass)=0. sec%transport(9,jclass)=0. ENDIF IF( zsum(2,jclass) .NE. 0 )THEN sec%transport( 4,jclass)=sec%transport( 4,jclass)+zsum( 4,jclass)/zsum(2,jclass) sec%transport( 6,jclass)=sec%transport( 6,jclass)+zsum( 6,jclass)/zsum(2,jclass) sec%transport( 8,jclass)=sec%transport( 8,jclass)+zsum( 8,jclass)/zsum(2,jclass) sec%transport(10,jclass)=sec%transport(10,jclass)+zsum(10,jclass)/zsum(2,jclass) ELSE sec%transport( 4,jclass)=0. sec%transport( 6,jclass)=0. sec%transport( 8,jclass)=0. sec%transport(10,jclass)=0. ENDIF ELSE sec%transport( 3,jclass)=0. sec%transport( 4,jclass)=0. sec%transport( 5,jclass)=0. sec%transport( 6,jclass)=0. sec%transport( 7,jclass)=0. sec%transport( 8,jclass)=0. sec%transport(10,jclass)=0. ENDIF ENDDO IF( sec%ll_ice_section ) THEN sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6 sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6 sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6 sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6 ENDIF END SUBROUTINE transport SUBROUTINE dia_dct_wri(kt,jsec,sec) !!------------------------------------------------------------- !! Write transport output in numdct !! !! Author: 10/05 Matthieu Laborie !!------------------------------------------------------------- !!arguments INTEGER, INTENT(IN) :: kt TYPE(SECTION), INTENT(INOUT) :: sec INTEGER ,INTENT(IN) :: jsec !!local declarations REAL(wp) ,DIMENSION(2):: zsumclass INTEGER :: jcl,ji !Dummy loop ! CHARACTER(len=35) :: zcritere CHARACTER(len=2) :: classe REAL(wp) :: zbnd1,zbnd2 !!------------------------------------------------------------- zsumclass(:)=0._wp DO jcl=1,MAX(1,sec%nb_class-1) ! Mean computation sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct) ! zcritere = 'no_classe' classe = 'N ' zbnd1 = 0._wp zbnd2 = 0._wp zsumclass(1)=zsumclass(1)+sec%transport(1,jcl) zsumclass(2)=zsumclass(2)+sec%transport(2,jcl) !insitu density classes transports IF( ( sec%zsigi(jcl) .NE. 99._wp ) .AND. & ( sec%zsigi(jcl+1) .NE. 99._wp ) )THEN classe = 'DI ' zbnd1 = sec%zsigi(jcl) zbnd2 = sec%zsigi(jcl+1) ENDIF !potential density classes transports IF( ( sec%zsigp(jcl) .NE. 99._wp ) .AND. & ( sec%zsigp(jcl+1) .NE. 99._wp ) )THEN classe = 'DP ' zbnd1 = sec%zsigp(jcl) zbnd2 = sec%zsigp(jcl+1) ENDIF !depth classes transports IF( ( sec%zlay(jcl) .NE. 99._wp ) .AND. & ( sec%zlay(jcl+1) .NE. 99._wp ) )THEN classe = 'Z ' zbnd1 = sec%zlay(jcl) zbnd2 = sec%zlay(jcl+1) ENDIF !salinity classes transports IF( ( sec%zsal(jcl) .NE. 99._wp ) .AND. & ( sec%zsal(jcl+1) .NE. 99._wp ) )THEN classe = 'S ' zbnd1 = sec%zsal(jcl) zbnd2 = sec%zsal(jcl+1) ENDIF !temperature classes transports IF( ( sec%ztem(jcl+1) .NE. 99. ) .AND. & ( sec%ztem(jcl+1) .NE. 99. ) ) THEN classe = 'T ' zbnd1 = sec%ztem(jcl) zbnd2 = sec%ztem(jcl+1) ENDIF ! ! Delate spaces and blank: ! zcritere = ADJUSTL(zcritere) ! ji = SCAN(TRIM(zcritere)," ") ! DO WHILE( ji .NE. 0 ) ! zcritere(ji:ji) = "_" ! ji = SCAN(TRIM(zcritere)," ") ! ENDDO !write volume transport per class WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,classe,zbnd1,zbnd2,& sec%transport(1,jcl),sec%transport(2,jcl), & sec%transport(1,jcl)+sec%transport(2,jcl) IF( sec%llstrpond )THEN !write heat transport per class WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,classe,zbnd1,zbnd2,& sec%transport(7,jcl)*rau0*rcp,sec%transport(8,jcl)* rau0*rcp, & ( sec%transport(7,jcl)+sec%transport(8,jcl) )*rau0*rcp !write salt transport per class WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,classe,zbnd1,zbnd2,& sec%transport(9,jcl)*rau0,sec%transport(10,jcl)*rau0,& (sec%transport(9,jcl)+sec%transport(10,jcl))*rau0 ENDIF ENDDO ! zcritere = "total" zbnd1 = 0._wp zbnd2 = 0._wp jcl=0 !write total volume transport WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,"total",zbnd1,zbnd2,& zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2) IF( sec%llstrpond )THEN !write total heat transport WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,"total",zbnd1,zbnd2,& zsumclass(7)* rau0*rcp,zsumclass(8)* rau0*rcp,(zsumclass(7)+zsumclass(8) )* rau0*rcp !write total salt transport WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & jcl,"total",zbnd1,zbnd2,& zsumclass(9)*rau0,zsumclass(10)*rau0,(zsumclass(9)+zsumclass(10))*rau0 ENDIF IF ( sec%ll_ice_section) THEN !write total ice volume transport WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,& jcl,"ice_vol",zbnd1,zbnd2,& sec%transport(9,1),sec%transport(10,1),& sec%transport(9,1)+sec%transport(10,1) !write total ice surface transport WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,& jcl,"ice_surf",zbnd1,zbnd2,& sec%transport(11,1),sec%transport(12,1), & sec%transport(11,1)+sec%transport(12,1) ENDIF !--------- !117 format means: X bound1 bound2 !with X=T/S/D/Z for Temperature/Salinity/Density/Depth classes !or X=N if there is No class !or X=total/ice_surf/ice_vol 118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4) 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) !--------- END SUBROUTINE dia_dct_wri FUNCTION interp(ki, kj, kk, cd_point, ptab) !!---------------------------------------------------------------------- !! !! !! !! !!---------------------------------------------------------------------- !*arguments INTEGER, INTENT(IN) :: ki, kj, kk CHARACTER(len=1), INTENT(IN) :: cd_point REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab REAL(wp) :: interp !*local declations INTEGER :: ii1, ij1, ii2, ij2 REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis !!---------------------------------------------------------------------- IF( cd_point=='U' )THEN ii1 = ki ; ij1 = kj ii2 = ki+1 ; ij2 = kj ELSE ! cd_point=='V' ii1 = ki ; ij1 = kj ii2 = ki ; ij2 = kj+1 ENDIF #if defined key_vvl ze3w = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk) zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) ) zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk) zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) ) zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk) #else ze3w = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk) zmax1 = ze3w / fse3w(ii2,ij2,kk) zmax2 = -ze3w / fse3w(ii1,ij1,kk) #endif IF(kk .NE. 1)THEN IF( ze3w >= 0. )THEN !zbis zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) ! result interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis ) ELSE !zbis zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) ! result interp = umask(ii1,ij1,kk) * 0.5*( zbis + ptab(ii2,ij2,kk) ) ENDIF ELSE interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) ) ENDIF END FUNCTION interp #else !!---------------------------------------------------------------------- !! Default option : Dummy module !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag CONTAINS SUBROUTINE dia_dct( kt ) ! Dummy routine INTEGER, INTENT( in ) :: kt ! ocean time-step index WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt END SUBROUTINE dia_dct #endif END MODULE diadct