- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/TOP/TRP/trcrad.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/TOP/TRP/trcrad.F90
r12178 r12928 6 6 !! History : - ! 01-01 (O. Aumont & E. Kestenare) Original code 7 7 !! 1.0 ! 04-03 (C. Ethe) free form F90 8 !! 4.1 ! 08-19 (A. Coward, D. Storkey) tidy up using new time-level indices 8 9 !!---------------------------------------------------------------------- 9 10 #if defined key_top … … 30 31 REAL(wp), DIMENSION(:,:), ALLOCATABLE:: gainmass 31 32 33 !! * Substitutions 34 # include "do_loop_substitute.h90" 32 35 !!---------------------------------------------------------------------- 33 36 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 37 40 CONTAINS 38 41 39 SUBROUTINE trc_rad( kt )42 SUBROUTINE trc_rad( kt, Kbb, Kmm, ptr ) 40 43 !!---------------------------------------------------------------------- 41 44 !! *** ROUTINE trc_rad *** … … 52 55 !! (the total CFC content is not strictly preserved) 53 56 !!---------------------------------------------------------------------- 54 INTEGER, INTENT(in) :: kt ! ocean time-step index 57 INTEGER, INTENT(in ) :: kt ! ocean time-step index 58 INTEGER, INTENT(in ) :: Kbb, Kmm ! time level indices 59 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation 55 60 ! 56 61 CHARACTER (len=22) :: charout … … 59 64 IF( ln_timing ) CALL timing_start('trc_rad') 60 65 ! 61 IF( ln_age ) CALL trc_rad_sms( kt, trb, trn, jp_age , jp_age ) ! AGE62 IF( ll_cfc ) CALL trc_rad_sms( kt, trb, trn, jp_cfc0, jp_cfc1 ) ! CFC model63 IF( ln_c14 ) CALL trc_rad_sms( kt, trb, trn, jp_c14 , jp_c14 ) ! C1464 IF( ln_pisces ) CALL trc_rad_sms( kt, trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model65 IF( ln_my_trc ) CALL trc_rad_sms( kt, trb, trn, jp_myt0, jp_myt1 ) ! MY_TRC model66 ! 67 IF( ln_ctl) THEN ! print mean trends (used for debugging)66 IF( ln_age ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_age , jp_age ) ! AGE 67 IF( ll_cfc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_cfc0, jp_cfc1 ) ! CFC model 68 IF( ln_c14 ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_c14 , jp_c14 ) ! C14 69 IF( ln_pisces ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model 70 IF( ln_my_trc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_myt0, jp_myt1 ) ! MY_TRC model 71 ! 72 IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging) 68 73 WRITE(charout, FMT="('rad')") 69 74 CALL prt_ctl_trc_info( charout ) 70 CALL prt_ctl_trc( tab4d= trn, mask=tmask, clinfo=ctrcnm )75 CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm ) 71 76 ENDIF 72 77 ! … … 87 92 !!---------------------------------------------------------------------- 88 93 ! 89 REWIND( numnat_ref ) ! namtrc_rad in reference namelist90 94 READ ( numnat_ref, namtrc_rad, IOSTAT = ios, ERR = 907) 91 95 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_rad in reference namelist' ) 92 REWIND( numnat_cfg ) ! namtrc_rad in configuration namelist93 96 READ ( numnat_cfg, namtrc_rad, IOSTAT = ios, ERR = 908 ) 94 97 908 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtrc_rad in configuration namelist' ) … … 113 116 114 117 115 SUBROUTINE trc_rad_sms( kt, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv ) 116 !!----------------------------------------------------------------------------- 117 !! *** ROUTINE trc_rad_sms *** 118 !! 119 !! ** Purpose : "crappy" routine to correct artificial negative 120 !! concentrations due to isopycnal scheme 121 !! 122 !! ** Method : 2 cases : 123 !! - Set negative concentrations to zero while computing 124 !! the corresponding tracer content that is added to the 125 !! tracers. Then, adjust the tracer concentration using 126 !! a multiplicative factor so that the total tracer 127 !! concentration is preserved. 128 !! - simply set to zero the negative CFC concentration 129 !! (the total content of concentration is not strictly preserved) 130 !!-------------------------------------------------------------------------------- 131 INTEGER , INTENT(in ) :: kt ! ocean time-step index 132 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 133 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) :: ptrb , ptrn ! before and now traceur concentration 134 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 135 ! 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 ) ) 146 ! 147 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 148 ! 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 161 ! 162 DO jk = 1, jpkm1 163 DO jj = 1, jpj 164 DO ji = 1, jpi 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 179 END DO 180 END DO 181 END DO 182 ! 183 IF( l_trdtrc ) THEN 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 273 ENDIF 274 ! 275 END DO 276 ! 277 ENDIF 118 SUBROUTINE trc_rad_sms( kt, Kbb, Kmm, ptr, jp_sms0, jp_sms1, cpreserv ) 119 !!----------------------------------------------------------------------------- 120 !! *** ROUTINE trc_rad_sms *** 121 !! 122 !! ** Purpose : "crappy" routine to correct artificial negative 123 !! concentrations due to isopycnal scheme 124 !! 125 !! ** Method : 2 cases : 126 !! - Set negative concentrations to zero while computing 127 !! the corresponding tracer content that is added to the 128 !! tracers. Then, adjust the tracer concentration using 129 !! a multiplicative factor so that the total tracer 130 !! concentration is preserved. 131 !! - simply set to zero the negative CFC concentration 132 !! (the total content of concentration is not strictly preserved) 133 !!-------------------------------------------------------------------------------- 134 INTEGER , INTENT(in ) :: kt ! ocean time-step index 135 INTEGER , INTENT(in ) :: Kbb, Kmm ! time level indices 136 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 137 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! before and now traceur concentration 138 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 139 ! 140 INTEGER :: ji, ji2, jj, jj2, jk, jn, jt ! dummy loop indices 141 INTEGER :: icnt, itime 142 LOGICAL :: lldebug = .FALSE. ! local logical 143 REAL(wp):: zcoef, zs2rdt, ztotmass 144 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 145 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 146 !!---------------------------------------------------------------------- 147 ! 148 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 149 zs2rdt = 1. / ( 2. * rn_Dt ) 150 ! 151 DO jt = 1,2 ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields 152 IF( jt == 1 ) itime = Kbb 153 IF( jt == 2 ) itime = Kmm 154 155 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 156 ! 157 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 158 159 DO jn = jp_sms0, jp_sms1 160 ztrneg(:,:,jn) = SUM( MIN( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 161 ztrpos(:,:,jn) = SUM( MAX( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 162 END DO 163 CALL sum3x3( ztrneg ) 164 CALL sum3x3( ztrpos ) 165 166 DO jn = jp_sms0, jp_sms1 167 ! 168 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr(:,:,:,:,Kbb) for trend computation 169 ! 170 DO_3D_11_11( 1, jpkm1 ) 171 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 172 ! 173 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * tmask(ji,jj,jk) ! really needed? 174 IF( ptr(ji,jj,jk,jn,itime) < 0. ) ptr(ji,jj,jk,jn,itime) = 0. ! suppress negative values 175 IF( ptr(ji,jj,jk,jn,itime) > 0. ) THEN ! use positive values to compensate mass gain 176 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptr > 0 177 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * zcoef 178 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 179 gainmass(jn,1) = gainmass(jn,1) - ptr(ji,jj,jk,jn,itime) * cvol(ji,jj,jk) ! we are adding mass... 180 ptr(ji,jj,jk,jn,itime) = 0. ! limit the compensation to keep positive value 181 ENDIF 182 ENDIF 183 ! 184 ENDIF 185 END_3D 186 ! 187 IF( l_trdtrc ) THEN 188 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 189 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 190 ENDIF 191 ! 192 END DO 193 194 IF( kt == nitend ) THEN 195 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 196 DO jn = jp_sms0, jp_sms1 197 IF( gainmass(jn,1) > 0. ) THEN 198 ztotmass = glob_sum( 'trcrad', ptr(:,:,:,jn,itime) * cvol(:,:,:) ) 199 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 200 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 201 END IF 202 END DO 203 ENDIF 204 205 DEALLOCATE( ztrneg, ztrpos ) 206 ! 207 ELSE !== total CFC content is NOT strictly preserved ==! 208 ! 209 DO jn = jp_sms0, jp_sms1 210 ! 211 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr for trend computation 212 ! 213 WHERE( ptr(:,:,:,jn,itime) < 0. ) ptr(:,:,:,jn,itime) = 0. 214 ! 215 IF( l_trdtrc ) THEN 216 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 217 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 218 ENDIF 219 ! 220 END DO 221 ! 222 ENDIF 223 ! 224 END DO 278 225 ! 279 226 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) … … 286 233 !!---------------------------------------------------------------------- 287 234 CONTAINS 288 SUBROUTINE trc_rad( kt ) ! Empty routine235 SUBROUTINE trc_rad( kt, Kbb, Kmm ) ! Empty routine 289 236 INTEGER, INTENT(in) :: kt 237 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 290 238 WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt 291 239 END SUBROUTINE trc_rad
Note: See TracChangeset
for help on using the changeset viewer.
