Changeset 786 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
- Timestamp:
- 2008-01-10T18:11:23+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
r719 r786 4 4 !! Ocean active tracers: horizontal component of the lateral tracer mixing trend 5 5 !!============================================================================== 6 !! History : 8.0 ! 97-07 (G. Madec) Original code 7 !! NEMO ! 02-08 (G. Madec) F90: Free form and module 8 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC 9 !!---------------------------------------------------------------------- 6 10 #if defined key_ldfslp || defined key_esopa 7 11 !!---------------------------------------------------------------------- … … 12 16 !! ldfght : ??? 13 17 !!---------------------------------------------------------------------- 14 !! * Modules used15 USE oce ! ocean dynamics and tracers variables16 18 USE dom_oce ! ocean space and time domain variables 17 19 USE ldftra_oce ! ocean active tracers: lateral physics … … 42 44 CONTAINS 43 45 44 SUBROUTINE tra_ldf_bilapg( kt ) 46 SUBROUTINE tra_ldf_bilapg( kt, cdtype, ktra, pgtu, pgtv, & 47 & ptb , pta ) 45 48 !!---------------------------------------------------------------------- 46 49 !! *** ROUTINE tra_ldf_bilapg *** … … 55 58 !! -1- compute the geopotential harmonic operator applied to 56 59 !! (tb,sb) and multiply it by the eddy diffusivity coefficient 57 !! (done by a call to ldfght routine, result in (wk1,wk2) arrays).60 !! (done by a call to ldfght routine, result in wk1 array). 58 61 !! Applied the domain lateral boundary conditions by call to lbc_lnk 59 62 !! -2- compute the geopotential harmonic operator applied to 60 !! (wk1,wk2) by a second call to ldfght routine (result in (wk3,wk4)63 !! wk1 by a second call to ldfght routine (result in wk2) 61 64 !! arrays). 62 !! -3- Add this trend to the general trend (ta,sa):63 !! (ta,sa) = (ta,sa) + (wk3,wk4)64 !! 65 !! ** Action : - Update (ta,sa)arrays with the before geopotential65 !! -3- Add this trend to the general trend pta: 66 !! pta = pta + wk2 67 !! 68 !! ** Action : - Update pta arrays with the before geopotential 66 69 !! biharmonic mixing trend. 67 70 !! … … 71 74 !! 9.0 ! 04-08 (C. Talandier) New trends organization 72 75 !!---------------------------------------------------------------------- 73 !! * Modules used 74 USE oce , wk1 => ua, & ! use ua as workspace 75 & wk2 => va ! use va as workspace 76 77 !! * Arguments 78 INTEGER, INTENT( in ) :: kt ! ocean time-step index 79 80 !! * Local declarations 76 INTEGER , INTENT(in ) :: kt ! ocean time-step index 77 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 78 INTEGER , INTENT(in ) :: ktra ! tracer index 79 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj) :: pgtu, pgtv ! tracer gradient at pstep levels 80 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb ! before tracer field 81 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 82 !! 81 83 INTEGER :: ji, jj, jk ! dummy loop indices 82 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 83 wk3, wk4 ! work array used for rotated biharmonic 84 ! ! operator on tracers and/or momentum 84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: wk1, wk2 ! workspace arrays 85 85 !!---------------------------------------------------------------------- 86 86 … … 91 91 ENDIF 92 92 93 ! 1. Laplacian of (tb,sb) * aht 94 ! ----------------------------- 95 ! rotated harmonic operator applied to (tb,sb) 96 ! and multiply by aht (output in (wk1,wk2) ) 97 98 CALL ldfght ( kt, tb, sb, wk1, wk2, 1 ) 99 100 101 ! Lateral boundary conditions on (wk1,wk2) (unchanged sign) 102 CALL lbc_lnk( wk1, 'T', 1. ) ; CALL lbc_lnk( wk2, 'T', 1. ) 103 104 ! 2. Bilaplacian of (tb,sb) 105 ! ------------------------- 106 ! rotated harmonic operator applied to (wk1,wk2) 107 ! (output in (wk3,wk4) ) 108 109 CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 ) 110 111 112 ! 3. Update the tracer trends (j-slab : 2, jpj-1) 113 ! --------------------------- 114 ! ! =============== 115 DO jj = 2, jpjm1 ! Vertical slab 116 ! ! =============== 117 DO jk = 1, jpkm1 93 ! Laplacian of ptb * aht 94 95 CALL ldfght ( kt, cdtype, ktra, ptb, wk1, 1 ) ! rotated laplacian applied to ptb and * aht (output in wk1 ) 96 97 CALL lbc_lnk( wk1, 'T', 1. ) ! Lateral boundary conditions on wk1 (unchanged sign) 98 99 ! Bilaplacian of ptb 100 101 CALL ldfght ( kt, cdtype, ktra, wk1, wk2, 2 ) ! rotated laplacian applied to wk1 (output in wk2 ) 102 103 104 ! Update the tracer trends 105 DO jk = 1, jpkm1 106 DO jj = 2, jpjm1 118 107 DO ji = 2, jpim1 119 ! add it to the general tracer trends 120 ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk) 121 sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk) 122 END DO 123 END DO 124 ! ! =============== 125 END DO ! End of slab 126 ! ! =============== 127 108 pta(ji,jj,jk) = pta(ji,jj,jk) + wk2(ji,jj,jk) 109 END DO 110 END DO 111 END DO 112 ! 128 113 END SUBROUTINE tra_ldf_bilapg 129 114 130 115 131 SUBROUTINE ldfght ( kt, pt, ps, plt, pls, kaht )116 SUBROUTINE ldfght ( kt, cdtype, ktra, pt, plt, kaht ) 132 117 !!---------------------------------------------------------------------- 133 118 !! *** ROUTINE ldfght *** 134 119 !! 135 !! ** Purpose : Apply a geopotential harmonic operator to (pt,ps)and120 !! ** Purpose : Apply a geopotential harmonic operator to p and 136 121 !! multiply it by the eddy diffusivity coefficient (if kaht=1). 137 122 !! Routine only used in s-coordinates (l_sco=T) with bilaplacian … … 140 125 !! 141 126 !! ** Method : The harmonic operator rotated along geopotential 142 !! surfaces is applied to (pt,ps)using the slopes of geopotential127 !! surfaces is applied to pt using the slopes of geopotential 143 128 !! surfaces computed in inildf routine. The result is provided in 144 !! (plt,pls)arrays. It is computed in 2 steps:129 !! plt arrays. It is computed in 2 steps: 145 130 !! 146 131 !! First step: horizontal part of the operator. It is computed on 147 !! ========== pt as follows (idem on ps)132 !! ========== pt as follows 148 133 !! horizontal fluxes : 149 134 !! zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ] … … 154 139 !! 155 140 !! Second step: vertical part of the operator. It is computed on 156 !! =========== pt as follows (idem on ps)141 !! =========== pt as follows 157 142 !! vertical fluxes : 158 143 !! zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2) dk-1[ pt ] … … 168 153 !! * Action : 169 154 !! 'key_trdtra' defined: the trend is saved for diagnostics. 170 !! 171 !! History : 172 !! 8.0 ! 97-07 (G. Madec) Original code 173 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 174 !!---------------------------------------------------------------------- 175 !! * Arguments 176 INTEGER, INTENT( in ) :: kt ! ocean time-step index 177 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & 178 pt, ps ! tracer fields (before t and s for 1st call 179 ! ! and laplacian of these fields for 2nd call. 180 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & 181 plt, pls ! partial harmonic operator applied to 182 ! ! pt & ps components except 183 ! ! second order vertical derivative term) 184 INTEGER, INTENT( in ) :: & 185 kaht ! =1 multiply the laplacian by the eddy diffusivity coeff. 186 ! ! =2 no multiplication 187 188 !! * Local declarations 189 INTEGER :: ji, jj, jk ! dummy loop indices 190 REAL(wp) :: & 191 zabe1, zabe2, zmku, zmkv, & ! temporary scalars 192 zbtr, ztah, zsah, ztav, zsav, & 193 zcof0, zcof1, zcof2, & 194 zcof3, zcof4 195 REAL(wp), DIMENSION(jpi,jpj) :: & 196 zftu, zfsu, & ! workspace 197 zdkt, zdk1t, & 198 zdks, zdk1s 199 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 200 zftv, zfsv ! workspace (only v components for ptr) 201 REAL(wp), DIMENSION(jpi,jpk) :: & 202 zftw, zfsw, & ! workspace 203 zdit, zdjt, zdj1t, & 204 zdis, zdjs, zdj1s 155 !!---------------------------------------------------------------------- 156 INTEGER , INTENT(in ) :: kt ! ocean time-step index 157 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 158 INTEGER , INTENT(in ) :: ktra ! tracer index 159 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pt ! before tracer field (1st call) 160 ! ! laplacian of the tracer field (2nd call) 161 REAL(wp) , INTENT(out), DIMENSION(jpi,jpj,jpk) :: plt ! harmonic operator applied to pt 162 INTEGER , INTENT(in ) :: kaht ! =1 multiply plt by aht 163 ! ! =2 no multiplication 164 !! 165 INTEGER :: ji, jj, jk ! dummy loop indices 166 REAL(wp) :: zabe1, zabe2, zmku, zmkv ! temporary scalars 167 REAL(wp) :: zbtr, ztah, ztav 168 REAL(wp) :: zcof0, zcof1, zcof2 169 REAL(wp) :: zcof3, zcof4 170 REAL(wp), DIMENSION(jpi,jpj) :: zftu, zdkt, zdk1t ! 2D workspace 171 REAL(wp), DIMENSION(jpi,jpk) :: zftw, zdit, zdjt, zdj1t ! 2D workspace 172 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zftv ! workspace (only v components for ptr) 205 173 !!---------------------------------------------------------------------- 206 174 … … 209 177 ! ! ********** ! ! =============== 210 178 211 ! I.1 Vertical gradient of pt a nd ps at level jk and jk+1212 ! ------------------------------------------------ -------179 ! I.1 Vertical gradient of pt at level jk and jk+1 180 ! ------------------------------------------------ 213 181 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 214 182 215 183 zdk1t(:,:) = ( pt(:,:,jk) - pt(:,:,jk+1) ) * tmask(:,:,jk+1) 216 zdk1s(:,:) = ( ps(:,:,jk) - ps(:,:,jk+1) ) * tmask(:,:,jk+1)217 184 218 185 IF( jk == 1 ) THEN 219 186 zdkt(:,:) = zdk1t(:,:) 220 zdks(:,:) = zdk1s(:,:)221 187 ELSE 222 188 zdkt(:,:) = ( pt(:,:,jk-1) - pt(:,:,jk) ) * tmask(:,:,jk) 223 zdks(:,:) = ( ps(:,:,jk-1) - ps(:,:,jk) ) * tmask(:,:,jk)224 189 ENDIF 225 190 … … 250 215 + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) & 251 216 +zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) 252 253 zfsu(ji,jj)= umask(ji,jj,jk) * &254 ( zabe1 *( ps(ji+1,jj,jk) - ps(ji,jj,jk) ) &255 + zcof1 *( zdks (ji+1,jj) + zdk1s(ji,jj) &256 +zdk1s(ji+1,jj) + zdks (ji,jj) ) )257 258 zfsv(ji,jj,jk)= vmask(ji,jj,jk) * &259 ( zabe2 *( ps(ji,jj+1,jk) - ps(ji,jj,jk) ) &260 + zcof2 *( zdks (ji,jj+1) + zdk1s(ji,jj) &261 +zdk1s(ji,jj+1) + zdks (ji,jj) ) )262 217 END DO 263 218 END DO … … 270 225 DO ji = 2 , jpim1 271 226 ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) 272 zsah = zfsu(ji,jj) - zfsu(ji-1,jj) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk)273 227 plt(ji,jj,jk) = ztah 274 pls(ji,jj,jk) = zsah275 228 END DO 276 229 END DO … … 279 232 ! ! =============== 280 233 281 !!but this should be done somewhere after 282 ! "zonal" mean diffusive heat and salt transport 283 IF( ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 284 pht_ldf(:) = ptr_vj( zftv(:,:,:) ) 285 pst_ldf(:) = ptr_vj( zfsv(:,:,:) ) 234 ! "Poleward" diffusive heat or salt transport 235 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 236 IF( ktra == jp_tem) pht_ldf(:) = ptr_vj( zftv(:,:,:) ) 237 IF( ktra == jp_sal) pst_ldf(:) = ptr_vj( zftv(:,:,:) ) 286 238 ENDIF 287 239 … … 296 248 DO ji = 1, jpim1 297 249 zdit (ji,jk) = ( pt(ji+1,jj ,jk) - pt(ji,jj ,jk) ) * umask(ji,jj ,jk) 298 zdis (ji,jk) = ( ps(ji+1,jj ,jk) - ps(ji,jj ,jk) ) * umask(ji,jj ,jk)299 250 zdjt (ji,jk) = ( pt(ji ,jj+1,jk) - pt(ji,jj ,jk) ) * vmask(ji,jj ,jk) 300 zdjs (ji,jk) = ( ps(ji ,jj+1,jk) - ps(ji,jj ,jk) ) * vmask(ji,jj ,jk)301 251 zdj1t(ji,jk) = ( pt(ji ,jj ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk) 302 zdj1s(ji,jk) = ( ps(ji ,jj ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)303 252 END DO 304 253 END DO … … 310 259 ! Surface and bottom vertical fluxes set to zero 311 260 zftw(:, 1 ) = 0.e0 312 zfsw(:, 1 ) = 0.e0313 261 zftw(:,jpk) = 0.e0 314 zfsw(:,jpk) = 0.e0315 262 316 263 ! interior (2=<jk=<jpk-1) … … 336 283 + zcof4 * ( zdjt (ji ,jk-1) + zdj1t(ji ,jk) & 337 284 +zdj1t(ji ,jk-1) + zdjt (ji ,jk) ) ) 338 339 zfsw(ji,jk) = tmask(ji,jj,jk) * &340 ( zcof0 * ( ps (ji,jj,jk-1) - ps (ji,jj,jk) ) &341 + zcof3 * ( zdis (ji ,jk-1) + zdis (ji-1,jk) &342 +zdis (ji-1,jk-1) + zdis (ji ,jk) ) &343 + zcof4 * ( zdjs (ji ,jk-1) + zdj1s(ji ,jk) &344 +zdj1s(ji ,jk-1) + zdjs (ji ,jk) ) )345 285 END DO 346 286 END DO … … 358 298 ! vertical divergence 359 299 ztav = zftw(ji,jk) - zftw(ji,jk+1) 360 zsav = zfsw(ji,jk) - zfsw(ji,jk+1)361 300 ! harmonic operator applied to (pt,ps) and multiply by aht 362 301 plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr 363 pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr364 302 END DO 365 303 END DO … … 372 310 ! vertical divergence 373 311 ztav = zftw(ji,jk) - zftw(ji,jk+1) 374 zsav = zfsw(ji,jk) - zfsw(ji,jk+1)375 312 ! harmonic operator applied to (pt,ps) 376 313 plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr 377 pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr378 314 END DO 379 315 END DO
Note: See TracChangeset
for help on using the changeset viewer.
