New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10425 for NEMO/trunk/src/TOP/TRP/trcrad.F90 – NEMO

Ignore:
Timestamp:
2018-12-19T22:54:16+01:00 (5 years ago)
Author:
smasson
Message:

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/TOP/TRP/trcrad.F90

    r10068 r10425  
    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( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 
    153             ztrcorn = glob_sum( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 
    154             ! 
    155             ztrmasb = glob_sum( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 
    156             ztrmasn = glob_sum( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 
    157             ! 
    158             IF( ztrcorb /= 0 ) THEN 
    159                zcoef = 1. + ztrcorb / ztrmasb 
    160                DO jk = 1, jpkm1 
    161                   ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) ) 
    162                   ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk) 
    163                END DO 
    164             ENDIF 
    165             ! 
    166             IF( ztrcorn /= 0 ) THEN 
    167                zcoef = 1. + ztrcorn / ztrmasn 
    168                DO jk = 1, jpkm1 
    169                   ptrn(:,:,jk,jn) = MAX( 0., ptrn(:,:,jk,jn) ) 
    170                   ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk) 
    171                END DO 
    172             ENDIF 
    173             ! 
    174             IF( l_trdtrc ) THEN 
    175                ! 
    176                zs2rdt = 1. / ( 2. * rdt ) 
    177                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    178                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
    179                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
    180                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    181               ! 
    182             ENDIF 
    183             ! 
    184          END DO 
    185          ! 
    186       ELSE                                !==  total CFC content is NOT strictly preserved  ==! 
    187          ! 
    188          DO jn = jp_sms0, jp_sms1   
    189             ! 
    190             IF( l_trdtrc ) THEN 
    191                ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation 
    192                ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation 
    193             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            
    194161            ! 
    195162            DO jk = 1, jpkm1 
    196163               DO jj = 1, jpj 
    197164                  DO ji = 1, jpi 
    198                      ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) ) 
    199                      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 
    200179                  END DO 
    201180               END DO 
    … …  
    203182            ! 
    204183            IF( l_trdtrc ) THEN 
    205                ! 
    206                zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 
    207                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    208                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
    209                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
    210                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    211               ! 
     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 
    212273            ENDIF 
    213274            ! 
    … …  
    216277      ENDIF 
    217278      ! 
    218       IF( l_trdtrc )  DEALLOCATE( ztrtrdb, ztrtrdn ) 
     279      IF( l_trdtrc )  DEALLOCATE( ztrtrd ) 
    219280      ! 
    220281   END SUBROUTINE trc_rad_sms 
Note: See TracChangeset for help on using the changeset viewer.