- Timestamp:
- 2020-12-02T16:13:45+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13559sette10 ^/utils/CI/sette_MPI3_LoopFusion@13943 sette
-
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/TRA/traldf_triad.F90
r13899 r14012 13 13 USE oce ! ocean dynamics and active tracers 14 14 USE dom_oce ! ocean space and time domain 15 ! TEMP: [tiling] This change not necessary if XIOS has subdomain support 16 USE domain, ONLY : dom_tile 17 USE domutl, ONLY : is_tile 15 18 USE phycst ! physical constants 16 19 USE trc_oce ! share passive tracers/Ocean variables … … 33 36 PUBLIC tra_ldf_triad ! routine called by traldf.F90 34 37 35 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels36 37 38 LOGICAL :: l_ptr ! flag to compute poleward transport 38 39 LOGICAL :: l_hst ! flag to compute heat transport … … 49 50 CONTAINS 50 51 51 SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, & 52 & pgu , pgv , pgui, pgvi , & 53 & pt , pt2, pt_rhs, kjpt, kpass ) 52 SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, & 53 & pgu , pgv , pgui, pgvi, & 54 & pt, pt2, pt_rhs, kjpt, kpass ) 55 !! 56 INTEGER , INTENT(in ) :: kt ! ocean time-step index 57 INTEGER , INTENT(in ) :: kit000 ! first time step index 58 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 59 INTEGER , INTENT(in ) :: kjpt ! number of tracers 60 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 61 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 62 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 63 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 64 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 65 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 66 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 67 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pt_rhs ! tracer trend 68 !! 69 CALL tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu), & 70 & pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), & 71 & pt, is_tile(pt), pt2, is_tile(pt2), pt_rhs, is_tile(pt_rhs), kjpt, kpass ) 72 END SUBROUTINE tra_ldf_triad 73 74 75 SUBROUTINE tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah, & 76 & pgu , pgv , ktg , pgui, pgvi, ktgi, & 77 & pt, ktt, pt2, ktt2, pt_rhs, ktt_rhs, kjpt, kpass ) 54 78 !!---------------------------------------------------------------------- 55 79 !! *** ROUTINE tra_ldf_triad *** … … 77 101 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 78 102 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 79 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 80 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 81 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 103 INTEGER , INTENT(in ) :: ktah, ktg, ktgi, ktt, ktt2, ktt_rhs 104 REAL(wp), DIMENSION(A2D_T(ktah), JPK) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 105 REAL(wp), DIMENSION(A2D_T(ktg), KJPT), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 106 REAL(wp), DIMENSION(A2D_T(ktgi), KJPT), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 107 REAL(wp), DIMENSION(A2D_T(ktt), JPK,KJPT), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 108 REAL(wp), DIMENSION(A2D_T(ktt2), JPK,KJPT), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 109 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 85 110 ! 86 111 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 94 119 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 95 120 REAL(wp) :: zah, zah_slp, zaei_slp 96 REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw ! 3D - 121 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 122 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 123 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk) :: zdit, zdjt, zftu, zftv, ztfw ! 3D - 124 ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support 125 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 98 126 !!---------------------------------------------------------------------- 99 127 ! 100 IF( .NOT.ALLOCATED(zdkt3d) ) THEN 101 ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr ) 102 CALL mpp_sum ( 'traldf_triad', ierr ) 103 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays') 104 ENDIF 105 ! 106 IF( kpass == 1 .AND. kt == kit000 ) THEN 107 IF(lwp) WRITE(numout,*) 108 IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype 109 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 110 ENDIF 111 ! 112 l_hst = .FALSE. 113 l_ptr = .FALSE. 114 IF( cdtype == 'TRA' ) THEN 115 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 116 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 128 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 129 IF( kpass == 1 .AND. kt == kit000 ) THEN 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype 132 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 133 ENDIF 134 ! 135 l_hst = .FALSE. 136 l_ptr = .FALSE. 137 IF( cdtype == 'TRA' ) THEN 138 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 139 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 140 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 141 ENDIF 118 142 ENDIF 119 143 ! … … 128 152 IF( kpass == 1 ) THEN !== first pass only and whatever the tracer is ==! 129 153 ! 130 akz (:,:,:) = 0._wp 131 ah_wslp2(:,:,:) = 0._wp 132 IF( ln_ldfeiv_dia ) THEN 133 zpsi_uw(:,:,:) = 0._wp 134 zpsi_vw(:,:,:) = 0._wp 135 ENDIF 154 DO_3D( 0, 0, 0, 0, 1, jpk ) 155 akz (ji,jj,jk) = 0._wp 156 ah_wslp2(ji,jj,jk) = 0._wp 157 END_3D 136 158 ! 137 159 DO ip = 0, 1 ! i-k triads 138 160 DO kp = 0, 1 139 DO_3D( 1, 0, 1, 0, 1, jpkm1 )140 ze3wr = 1._wp / e3w(ji +ip,jj,jk+kp,Kmm)141 zbu = e1e2u(ji ,jj) * e3u(ji,jj,jk,Kmm)142 zah = 0.25_wp * pahu(ji ,jj,jk)143 zslope_skew = triadi_g(ji +ip,jj,jk,1-ip,kp)161 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 162 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 163 zbu = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm) 164 zah = 0.25_wp * pahu(ji-ip,jj,jk) 165 zslope_skew = triadi_g(ji,jj,jk,1-ip,kp) 144 166 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 145 zslope2 = zslope_skew + ( gdept(ji +1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)167 zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 146 168 zslope2 = zslope2 *zslope2 147 ah_wslp2(ji +ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2148 akz (ji +ip,jj,jk+kp) = akz (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj) &149 & * r1_e1u(ji ,jj) * umask(ji,jj,jk+kp)169 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 170 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj) & 171 & * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 150 172 ! 151 IF( ln_ldfeiv_dia ) zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &152 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew153 173 END_3D 154 174 END DO … … 157 177 DO jp = 0, 1 ! j-k triads 158 178 DO kp = 0, 1 159 DO_3D( 1, 0, 1, 0, 1, jpkm1 )160 ze3wr = 1.0_wp / e3w(ji,jj +jp,jk+kp,Kmm)161 zbv = e1e2v(ji,jj ) * e3v(ji,jj,jk,Kmm)162 zah = 0.25_wp * pahv(ji,jj ,jk)163 zslope_skew = triadj_g(ji,jj +jp,jk,1-jp,kp)179 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 180 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 181 zbv = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm) 182 zah = 0.25_wp * pahv(ji,jj-jp,jk) 183 zslope_skew = triadj_g(ji,jj,jk,1-jp,kp) 164 184 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 165 185 ! (do this by *adding* gradient of depth) 166 zslope2 = zslope_skew + ( gdept(ji,jj +1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)186 zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 167 187 zslope2 = zslope2 * zslope2 168 ah_wslp2(ji,jj +jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2169 akz (ji,jj +jp,jk+kp) = akz (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj) &170 & * r1_e2v(ji,jj ) * vmask(ji,jj,jk+kp)188 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 189 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp) & 190 & * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 171 191 ! 172 IF( ln_ldfeiv_dia ) zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &173 & + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew174 192 END_3D 175 193 END DO … … 179 197 ! 180 198 IF( ln_traldf_blp ) THEN ! bilaplacian operator 181 DO_3D( 1, 0, 1, 0, 2, jpkm1 )199 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 182 200 akz(ji,jj,jk) = 16._wp & 183 201 & * ah_wslp2 (ji,jj,jk) & … … 187 205 END_3D 188 206 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 189 DO_3D( 1, 0, 1, 0, 2, jpkm1 )207 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 190 208 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 191 209 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 195 213 ! 196 214 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 197 akz(:,:,:) = ah_wslp2(:,:,:) 198 ENDIF 199 ! 200 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 215 DO_3D( 0, 0, 0, 0, 1, jpk ) 216 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 217 END_3D 218 ENDIF 219 ! 220 ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 221 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 222 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 223 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 224 225 zpsi_uw(:,:,:) = 0._wp 226 zpsi_vw(:,:,:) = 0._wp 227 228 DO jp = 0, 1 229 DO kp = 0, 1 230 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 231 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 232 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 233 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 234 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 235 END_3D 236 END DO 237 END DO 238 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 239 240 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile ) 241 ENDIF 242 ENDIF 201 243 ! 202 244 ENDIF !== end 1st pass only ==! … … 234 276 DO jk = 1, jpkm1 235 277 ! !== Vertical tracer gradient at level jk and jk+1 236 zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 278 DO_2D( 1, 1, 1, 1 ) 279 zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 280 END_2D 237 281 ! 238 282 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 239 283 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 240 ELSE ; zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk) 284 ELSE 285 DO_2D( 1, 1, 1, 1 ) 286 zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 287 END_2D 241 288 ENDIF 242 289 ! … … 380 427 END DO ! end tracer loop 381 428 ! ! =============== 382 END SUBROUTINE tra_ldf_triad 429 END SUBROUTINE tra_ldf_triad_t 383 430 384 431 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.