Changeset 10408


Ignore:
Timestamp:
2018-12-18T12:21:21+01:00 (20 months ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 8 introduce sum3x3 and supress glob_sum in trcrad, see #2133

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90

    r10397 r10408  
    2121   USE in_out_manager  ! I/O manager 
    2222   USE lib_mpp         ! distributed memory computing 
     23   USE lbclnk          ! ocean lateral boundary conditions 
    2324 
    2425   IMPLICIT NONE 
     
    2829   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos) 
    2930   PUBLIC   local_sum     ! used in trcrad, local operation before glob_sum_delay 
     31   PUBLIC   sum3x3        ! used in trcrad, do a sum over 3x3 boxes 
    3032   PUBLIC   DDPDD         ! also used in closea module 
    3133   PUBLIC   glob_min, glob_max 
     
    4244   INTERFACE local_sum 
    4345      MODULE PROCEDURE local_sum_2d, local_sum_3d 
     46   END INTERFACE 
     47   INTERFACE sum3x3 
     48      MODULE PROCEDURE sum3x3_2d, sum3x3_3d 
    4449   END INTERFACE 
    4550   INTERFACE glob_min 
     
    194199   END FUNCTION local_sum_3d 
    195200 
     201!                          ! FUNCTION sum3x3 ! 
     202 
     203   SUBROUTINE sum3x3_2d( p2d ) 
     204      !!----------------------------------------------------------------------- 
     205      !!                  ***  routine sum3x3_2d  *** 
     206      !! 
     207      !! ** Purpose : sum over 3x3 boxes 
     208      !!---------------------------------------------------------------------- 
     209      REAL(wp), DIMENSION (:,:), INTENT(inout) ::   p2d 
     210      ! 
     211      INTEGER ::   ji, ji2, jj, jj2     ! dummy loop indices 
     212      !!---------------------------------------------------------------------- 
     213      ! 
     214      IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' )  
     215      IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' )  
     216      ! 
     217      DO jj = 1, jpj 
     218         DO ji = 1, jpi  
     219            IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN   ! bottom left corber of a 3x3 box 
     220               ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box 
     221               jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box 
     222               IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain 
     223                  p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2)) 
     224               ENDIF 
     225            ENDIF 
     226         END DO 
     227      END DO 
     228      CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. ) 
     229      IF( nbondi /= -1 ) THEN 
     230         IF( MOD(mig(    1), 3) == 1 )   p2d(    1,:) = p2d(    2,:) 
     231         IF( MOD(mig(    1), 3) == 2 )   p2d(    2,:) = p2d(    1,:) 
     232      ENDIF 
     233      IF( nbondi /=  1 ) THEN 
     234         IF( MOD(mig(jpi-2), 3) == 1 )   p2d(  jpi,:) = p2d(jpi-1,:) 
     235         IF( MOD(mig(jpi-2), 3) == 0 )   p2d(jpi-1,:) = p2d(  jpi,:) 
     236      ENDIF 
     237      IF( nbondj /= -1 ) THEN 
     238         IF( MOD(mjg(    1), 3) == 1 )   p2d(:,    1) = p2d(:,    2) 
     239         IF( MOD(mjg(    1), 3) == 2 )   p2d(:,    2) = p2d(:,    1) 
     240      ENDIF 
     241      IF( nbondj /=  1 ) THEN 
     242         IF( MOD(mjg(jpj-2), 3) == 1 )   p2d(:,  jpj) = p2d(:,jpj-1) 
     243         IF( MOD(mjg(jpj-2), 3) == 0 )   p2d(:,jpj-1) = p2d(:,  jpj) 
     244      ENDIF 
     245      CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. ) 
     246 
     247   END SUBROUTINE sum3x3_2d 
     248 
     249   SUBROUTINE sum3x3_3d( p3d ) 
     250      !!----------------------------------------------------------------------- 
     251      !!                  ***  routine sum3x3_3d  *** 
     252      !! 
     253      !! ** Purpose : sum over 3x3 boxes 
     254      !!---------------------------------------------------------------------- 
     255      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   p3d 
     256      ! 
     257      INTEGER ::   ji, ji2, jj, jj2, jn     ! dummy loop indices 
     258      INTEGER ::   ipn                      ! Third dimension size 
     259      !!---------------------------------------------------------------------- 
     260      ! 
     261      IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' )  
     262      IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' )  
     263      ipn = SIZE(p3d,3) 
     264      ! 
     265      DO jn = 1, ipn 
     266         DO jj = 1, jpj 
     267            DO ji = 1, jpi  
     268               IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN   ! bottom left corber of a 3x3 box 
     269                  ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box 
     270                  jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box 
     271                  IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain 
     272                     p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn)) 
     273                  ENDIF 
     274               ENDIF 
     275            END DO 
     276         END DO 
     277      END DO 
     278      CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. ) 
     279      IF( nbondi /= -1 ) THEN 
     280         IF( MOD(mig(    1), 3) == 1 )   p3d(    1,:,:) = p3d(    2,:,:) 
     281         IF( MOD(mig(    1), 3) == 2 )   p3d(    2,:,:) = p3d(    1,:,:) 
     282      ENDIF 
     283      IF( nbondi /=  1 ) THEN 
     284         IF( MOD(mig(jpi-2), 3) == 1 )   p3d(  jpi,:,:) = p3d(jpi-1,:,:) 
     285         IF( MOD(mig(jpi-2), 3) == 0 )   p3d(jpi-1,:,:) = p3d(  jpi,:,:) 
     286      ENDIF 
     287      IF( nbondj /= -1 ) THEN 
     288         IF( MOD(mjg(    1), 3) == 1 )   p3d(:,    1,:) = p3d(:,    2,:) 
     289         IF( MOD(mjg(    1), 3) == 2 )   p3d(:,    2,:) = p3d(:,    1,:) 
     290      ENDIF 
     291      IF( nbondj /=  1 ) THEN 
     292         IF( MOD(mjg(jpj-2), 3) == 1 )   p3d(:,  jpj,:) = p3d(:,jpj-1,:) 
     293         IF( MOD(mjg(jpj-2), 3) == 0 )   p3d(:,jpj-1,:) = p3d(:,  jpj,:) 
     294      ENDIF 
     295      CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. ) 
     296 
     297   END SUBROUTINE sum3x3_3d 
     298 
    196299 
    197300   SUBROUTINE DDPDD( ydda, yddb ) 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcrad.F90

    r10365 r10408  
    1919   USE trdtra 
    2020   USE prtctl_trc          ! Print control for debbuging 
     21   USE lib_fortran 
    2122 
    2223   IMPLICIT NONE 
     
    2728 
    2829   LOGICAL , PUBLIC ::   ln_trcrad           !: flag to artificially correct negative concentrations 
     30   REAL(wp), DIMENSION(:,:), ALLOCATABLE::   gainmass 
    2931 
    3032   !!---------------------------------------------------------------------- 
     
    104106         ENDIF 
    105107      ENDIF 
     108      ! 
     109      ALLOCATE( gainmass(jptra,2) ) 
     110      gainmass(:,:) = 0. 
    106111      ! 
    107112   END SUBROUTINE trc_rad_ini 
     
    129134      CHARACTER( len = 1), OPTIONAL          , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not 
    130135      ! 
    131       INTEGER ::   ji, jj, jk, jn     ! dummy loop indices 
    132       LOGICAL ::   lldebug = .FALSE.           ! local logical 
    133       REAL(wp)::   ztrcorb, ztrmasb, zs2rdt    ! temporary scalars 
    134       REAL(wp)::   zcoef  , ztrcorn, ztrmasn   !    -         - 
    135       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays 
    136       !!---------------------------------------------------------------------- 
    137       ! 
    138       IF( l_trdtrc )   ALLOCATE( ztrtrdb(jpi,jpj,jpk), ztrtrdn(jpi,jpj,jpk) ) 
     136      INTEGER ::   ji, ji2, jj, jj2, jk, jn     ! dummy loop indices 
     137      INTEGER ::   icnt 
     138      LOGICAL ::   lldebug = .FALSE.            ! local logical 
     139      REAL(wp)::   zcoef, zs2rdt, ztotmass 
     140      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrneg, ztrpos 
     141      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd   ! workspace arrays 
     142      !!---------------------------------------------------------------------- 
     143      ! 
     144      IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 
     145      zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 
    139146      ! 
    140147      IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==! 
    141148         ! 
    142          DO jn = jp_sms0, jp_sms1 
    143             ! 
    144             ztrcorb = 0._wp   ;   ztrmasb = 0._wp 
    145             ztrcorn = 0._wp   ;   ztrmasn = 0._wp 
    146             ! 
    147             IF( l_trdtrc ) THEN 
    148                ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation 
    149                ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation 
    150             ENDIF 
    151             !                                                         ! sum over the global domain  
    152             ztrcorb = glob_sum( 'trcrad', MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 
    153             ztrcorn = glob_sum( 'trcrad', MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 
    154             ! 
    155             IF( ztrcorb /= 0 ) THEN 
    156                ztrmasb = glob_sum( 'trcrad', MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 
    157                zcoef = 1. + ztrcorb / ztrmasb 
    158                DO jk = 1, jpkm1 
    159                   ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) ) 
    160                   ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk) 
    161                END DO 
    162             ENDIF 
    163             ! 
    164             IF( ztrcorn /= 0 ) THEN 
    165                ztrmasn = glob_sum( 'trcrad', MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 
    166                zcoef = 1. + ztrcorn / ztrmasn 
    167                DO jk = 1, jpkm1 
    168                   ptrn(:,:,jk,jn) = MAX( 0., ptrn(:,:,jk,jn) ) 
    169                   ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk) 
    170                END DO 
    171             ENDIF 
    172             ! 
    173             IF( l_trdtrc ) THEN 
    174                ! 
    175                zs2rdt = 1. / ( 2. * rdt ) 
    176                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    177                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
    178                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
    179                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    180               ! 
    181             ENDIF 
    182             ! 
    183          END DO 
    184          ! 
    185       ELSE                                !==  total CFC content is NOT strictly preserved  ==! 
    186          ! 
    187          DO jn = jp_sms0, jp_sms1   
    188             ! 
    189             IF( l_trdtrc ) THEN 
    190                ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation 
    191                ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation 
    192             ENDIF 
     149         ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 
     150 
     151         DO jn = jp_sms0, jp_sms1 
     152            ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values 
     153            ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values 
     154         END DO 
     155         CALL sum3x3( ztrneg ) 
     156         CALL sum3x3( ztrpos ) 
     157          
     158         DO jn = jp_sms0, jp_sms1 
     159            ! 
     160            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input trb for trend computation            
    193161            ! 
    194162            DO jk = 1, jpkm1 
    195163               DO jj = 1, jpj 
    196164                  DO ji = 1, jpi 
    197                      ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) ) 
    198                      ptrb(ji,jj,jk,jn) = MAX( 0. , ptrb(ji,jj,jk,jn) ) 
     165                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box 
     166                        ! 
     167                        ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed? 
     168                        IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0.       ! supress negative values 
     169                        IF( ptrb(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain 
     170                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0 
     171                           ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef 
     172                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value 
     173                              gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass... 
     174                              ptrb(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value 
     175                           ENDIF 
     176                        ENDIF 
     177                        ! 
     178                     ENDIF 
    199179                  END DO 
    200180               END DO 
     
    202182            ! 
    203183            IF( l_trdtrc ) THEN 
    204                ! 
    205                zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 
    206                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    207                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
    208                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
    209                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    210               ! 
     184               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 
     185               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling 
     186            ENDIF 
     187            ! 
     188         END DO 
     189  
     190         IF( kt == nitend ) THEN 
     191            CALL mpp_sum( 'trcrad', gainmass(:,1) ) 
     192            DO jn = jp_sms0, jp_sms1 
     193               IF( gainmass(jn,1) > 0. ) THEN 
     194                  ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) ) 
     195                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  & 
     196                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1) 
     197               END IF 
     198            END DO 
     199         ENDIF 
     200 
     201         DO jn = jp_sms0, jp_sms1 
     202            ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values 
     203            ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values 
     204         END DO 
     205         CALL sum3x3( ztrneg ) 
     206         CALL sum3x3( ztrpos ) 
     207          
     208         DO jn = jp_sms0, jp_sms1 
     209            ! 
     210            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input trb for trend computation 
     211            ! 
     212            DO jk = 1, jpkm1 
     213               DO jj = 1, jpj 
     214                  DO ji = 1, jpi 
     215                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box 
     216                        ! 
     217                        ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed? 
     218                        IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0.       ! supress negative values 
     219                        IF( ptrn(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain 
     220                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0 
     221                           ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef 
     222                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value 
     223                              gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass... 
     224                              ptrn(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value 
     225                           ENDIF 
     226                        ENDIF 
     227                        ! 
     228                     ENDIF 
     229                  END DO 
     230               END DO 
     231            END DO 
     232            ! 
     233            IF( l_trdtrc ) THEN 
     234               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 
     235               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling 
     236            ENDIF 
     237            ! 
     238         END DO 
     239  
     240         IF( kt == nitend ) THEN 
     241            CALL mpp_sum( 'trcrad', gainmass(:,2) ) 
     242            DO jn = jp_sms0, jp_sms1 
     243               IF( gainmass(jn,2) > 0. ) THEN 
     244                  ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) ) 
     245                  WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  & 
     246                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1) 
     247               END IF 
     248            END DO 
     249         ENDIF 
     250 
     251         DEALLOCATE( ztrneg, ztrpos ) 
     252         ! 
     253      ELSE                                !==  total CFC content is NOT strictly preserved  ==! 
     254         ! 
     255         DO jn = jp_sms0, jp_sms1   
     256            ! 
     257            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation 
     258            ! 
     259            WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0. 
     260            ! 
     261            IF( l_trdtrc ) THEN 
     262               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 
     263               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling 
     264            ENDIF 
     265            ! 
     266            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation 
     267            ! 
     268            WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0. 
     269            ! 
     270            IF( l_trdtrc ) THEN 
     271               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 
     272               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling 
    211273            ENDIF 
    212274            ! 
     
    215277      ENDIF 
    216278      ! 
    217       IF( l_trdtrc )  DEALLOCATE( ztrtrdb, ztrtrdn ) 
     279      IF( l_trdtrc )  DEALLOCATE( ztrtrd ) 
    218280      ! 
    219281   END SUBROUTINE trc_rad_sms 
Note: See TracChangeset for help on using the changeset viewer.