[3611] | 1 | MODULE traldf_lap_tam |
---|
| 2 | #if defined key_tam |
---|
| 3 | !!============================================================================== |
---|
| 4 | !! *** MODULE traldf_lap *** |
---|
| 5 | !! Ocean active tracers: horizontal component of the lateral tracer mixing trend |
---|
| 6 | !! Tangent and adjoint module |
---|
| 7 | !!============================================================================== |
---|
| 8 | !! History : 9.0 ! ????? |
---|
| 9 | !! History of the T&A module: |
---|
| 10 | !! 9.0 ! 09-03 (F. Vigilant) original version |
---|
| 11 | !! NEMO 3.4 ! 12-07 (P.-A. Bouttier) |
---|
| 12 | !!---------------------------------------------------------------------- |
---|
| 13 | |
---|
| 14 | !!---------------------------------------------------------------------- |
---|
| 15 | !! tra_ldf_lap : update the tracer trend with the horizontal diffusion |
---|
| 16 | !! using a iso-level harmonic (laplacien) operator. |
---|
| 17 | !!---------------------------------------------------------------------- |
---|
| 18 | !! * Modules used |
---|
| 19 | USE par_oce |
---|
| 20 | USE oce_tam |
---|
| 21 | USE dom_oce |
---|
| 22 | USE ldftra_oce |
---|
| 23 | USE in_out_manager |
---|
| 24 | USE gridrandom |
---|
| 25 | USE dotprodfld |
---|
| 26 | USE tstool_tam |
---|
| 27 | USE paresp |
---|
| 28 | USE trc_oce |
---|
| 29 | USE lib_mpp |
---|
| 30 | USE timing |
---|
| 31 | USE wrk_nemo |
---|
| 32 | |
---|
| 33 | IMPLICIT NONE |
---|
| 34 | PRIVATE |
---|
| 35 | |
---|
| 36 | !! * Routine accessibility |
---|
| 37 | PUBLIC tra_ldf_lap_tan ! routine called by tradldf_tam.F90 |
---|
| 38 | PUBLIC tra_ldf_lap_adj ! routine called by tradldf_tam.F90 |
---|
| 39 | PUBLIC tra_ldf_lap_adj_tst ! routine called by tradldf_tam.F90 |
---|
| 40 | |
---|
| 41 | REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: e1ur, e2vr ! scale factor coefficients |
---|
| 42 | |
---|
| 43 | !! * Substitutions |
---|
| 44 | # include "domzgr_substitute.h90" |
---|
| 45 | # include "ldftra_substitute.h90" |
---|
| 46 | # include "vectopt_loop_substitute.h90" |
---|
| 47 | !!---------------------------------------------------------------------- |
---|
| 48 | !! OPA 9.0 , LOCEAN-IPSL (2005) |
---|
| 49 | !! $Id: traldf_lap.F90 1152 2008-06-26 14:11:13Z rblod $ |
---|
| 50 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
---|
| 51 | !!---------------------------------------------------------------------- |
---|
| 52 | |
---|
| 53 | CONTAINS |
---|
| 54 | |
---|
| 55 | SUBROUTINE tra_ldf_lap_tan( kt, kit000, cdtype, pgu_tl, pgv_tl, & |
---|
| 56 | & ptb_tl, pta_tl, kjpt ) |
---|
| 57 | !!---------------------------------------------------------------------- |
---|
| 58 | !! *** ROUTINE tra_ldf_lap_tan *** |
---|
| 59 | !! |
---|
| 60 | !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive |
---|
| 61 | !! trend and add it to the general trend of tracer equation. |
---|
| 62 | !! |
---|
| 63 | !! ** Method : Second order diffusive operator evaluated using before |
---|
| 64 | !! fields (forward time scheme). The horizontal diffusive trends of |
---|
| 65 | !! temperature (idem for salinity) is given by: |
---|
| 66 | !! difft = 1/(e1t*e2t*e3t) { di-1[ aht e2u*e3u/e1u di(tb) ] |
---|
| 67 | !! + dj-1[ aht e1v*e3v/e2v dj(tb) ] } |
---|
| 68 | !! Note: key_zco defined, the e3t=e3u=e3v, the trend becomes: |
---|
| 69 | !! difft = 1/(e1t*e2t) { di-1[ aht e2u/e1u di(tb) ] |
---|
| 70 | !! + dj-1[ aht e1v/e2v dj(tb) ] } |
---|
| 71 | !! Add this trend to the general tracer trend (ta,sa): |
---|
| 72 | !! (ta,sa) = (ta,sa) + ( difft , diffs ) |
---|
| 73 | !! |
---|
| 74 | !! ** Action : - Update (ta,sa) arrays with the before iso-level |
---|
| 75 | !! harmonic mixing trend. |
---|
| 76 | !! |
---|
| 77 | !! History : |
---|
| 78 | !! 1.0 ! 87-06 (P. Andrich, D. L Hostis) Original code |
---|
| 79 | !! ! 91-11 (G. Madec) |
---|
| 80 | !! ! 95-11 (G. Madec) suppress volumetric scale factors |
---|
| 81 | !! ! 96-01 (G. Madec) statement function for e3 |
---|
| 82 | !! 8.5 ! 02-06 (G. Madec) F90: Free form and module |
---|
| 83 | !! 9.0 ! 04-08 (C. Talandier) New trends organization |
---|
| 84 | !! ! 05-11 (G. Madec) add zps case |
---|
| 85 | !! History of the tangent routine |
---|
| 86 | !! 9.0 ! 03-09 (F. Vigilant) tangent of 9.0 |
---|
| 87 | !!---------------------------------------------------------------------- |
---|
| 88 | |
---|
| 89 | !! * Arguments |
---|
| 90 | INTEGER, INTENT( in ) :: kt ! ocean time-step index |
---|
| 91 | INTEGER , INTENT(in ) :: kit000 ! first time step index |
---|
| 92 | CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) |
---|
| 93 | INTEGER , INTENT(in ) :: kjpt ! number of tracers |
---|
| 94 | REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu_tl, pgv_tl ! tracer gradient at pstep levels |
---|
| 95 | REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb_tl ! before and now tracer fields |
---|
| 96 | REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta_tl ! tracer trend |
---|
| 97 | !! * Local declarations |
---|
| 98 | INTEGER :: ji, jj, jk, jn, ierr ! dummy loop indices |
---|
| 99 | INTEGER :: iku, ikv ! temporary integers |
---|
| 100 | REAL(wp) :: zabe1, zabe2, zbtr ! temporary scalars |
---|
| 101 | REAL(wp) :: ztatl, zsatl ! temporary scalars |
---|
| 102 | REAL(wp), POINTER, DIMENSION(:,:,:) :: ztutl, ztvtl ! 3D workspace |
---|
| 103 | !!---------------------------------------------------------------------- |
---|
| 104 | ! |
---|
| 105 | CALL timing_start('tra_ldf_lap_tan') |
---|
| 106 | ! |
---|
| 107 | CALL wrk_alloc(jpi,jpj,jpk, ztutl, ztvtl) |
---|
| 108 | ! |
---|
| 109 | IF( kt == nit000 ) THEN |
---|
| 110 | IF(lwp) WRITE(numout,*) |
---|
| 111 | IF(lwp) WRITE(numout,*) 'tra_ldf_lap_tan : iso-level laplacian diffusion on ', cdtype |
---|
| 112 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' |
---|
| 113 | IF( .NOT. ALLOCATED( e1ur ) ) THEN |
---|
| 114 | ! This routine may be called for both active and passive tracers. |
---|
| 115 | ! Allocate and set saved arrays on first call only. |
---|
| 116 | ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr ) |
---|
| 117 | IF( lk_mpp ) CALL mpp_sum( ierr ) |
---|
| 118 | IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' ) |
---|
| 119 | ! |
---|
| 120 | e1ur(:,:) = e2u(:,:) / e1u(:,:) |
---|
| 121 | e2vr(:,:) = e1v(:,:) / e2v(:,:) |
---|
| 122 | ENDIF |
---|
| 123 | ENDIF |
---|
| 124 | |
---|
| 125 | ! |
---|
| 126 | DO jn = 1, kjpt |
---|
| 127 | ! ! ============= |
---|
| 128 | DO jk = 1, jpkm1 ! Vertical slab |
---|
| 129 | ! ! ============= |
---|
| 130 | ! 1. First derivative (gradient) |
---|
| 131 | ! ------------------- |
---|
| 132 | DO jj = 1, jpjm1 |
---|
| 133 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
| 134 | zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk) |
---|
| 135 | zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk) |
---|
| 136 | ztutl(ji,jj,jk) = zabe1 * ( ptb_tl(ji+1,jj ,jk,jn) - ptb_tl(ji,jj,jk,jn) ) |
---|
| 137 | ztvtl(ji,jj,jk) = zabe2 * ( ptb_tl(ji ,jj+1,jk,jn) - ptb_tl(ji,jj,jk,jn) ) |
---|
| 138 | END DO |
---|
| 139 | END DO |
---|
| 140 | IF( ln_zps ) THEN ! set gradient at partial step level |
---|
| 141 | DO jj = 1, jpjm1 |
---|
| 142 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
| 143 | ! last level |
---|
| 144 | iku = mbku(ji,jj) |
---|
| 145 | ikv = mbku(ji,jj) |
---|
| 146 | IF( iku == jk ) THEN |
---|
| 147 | zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku) |
---|
| 148 | ztutl(ji,jj,jk) = zabe1 * pgu_tl(ji,jj,jn) |
---|
| 149 | ENDIF |
---|
| 150 | IF( ikv == jk ) THEN |
---|
| 151 | zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv) |
---|
| 152 | ztvtl(ji,jj,jk) = zabe2 * pgv_tl(ji,jj,jn) |
---|
| 153 | ENDIF |
---|
| 154 | END DO |
---|
| 155 | END DO |
---|
| 156 | ENDIF |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | ! 2. Second derivative (divergence) |
---|
| 160 | ! -------------------- |
---|
| 161 | DO jj = 2, jpjm1 |
---|
| 162 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
| 163 | zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)) |
---|
| 164 | ! horizontal diffusive trends |
---|
| 165 | ztatl = zbtr * ( ztutl(ji,jj,jk) - ztutl(ji-1,jj ,jk) & |
---|
| 166 | & + ztvtl(ji,jj,jk) - ztvtl(ji ,jj-1,jk) ) |
---|
| 167 | ! add it to the general tracer trends |
---|
| 168 | pta_tl(ji,jj,jk,jn) = pta_tl(ji,jj,jk,jn) + ztatl |
---|
| 169 | END DO |
---|
| 170 | END DO |
---|
| 171 | ! ! ============= |
---|
| 172 | END DO ! End of slab |
---|
| 173 | ! ! ============= |
---|
| 174 | END DO |
---|
| 175 | ! |
---|
| 176 | CALL timing_stop('tra_ldf_lap_tan') |
---|
| 177 | ! |
---|
| 178 | CALL wrk_dealloc(jpi,jpj,jpk, ztutl, ztvtl) |
---|
| 179 | ! |
---|
| 180 | |
---|
| 181 | END SUBROUTINE tra_ldf_lap_tan |
---|
| 182 | |
---|
| 183 | |
---|
| 184 | SUBROUTINE tra_ldf_lap_adj( kt, kit000, cdtype, pgu_ad, pgv_ad, & |
---|
| 185 | & ptb_ad, pta_ad, kjpt ) |
---|
| 186 | !!---------------------------------------------------------------------- |
---|
| 187 | !! *** ROUTINE tra_ldf_lap_adj *** |
---|
| 188 | !! |
---|
| 189 | !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive |
---|
| 190 | !! trend and add it to the general trend of tracer equation. |
---|
| 191 | !! |
---|
| 192 | !! ** Method : Second order diffusive operator evaluated using before |
---|
| 193 | !! fields (forward time scheme). The horizontal diffusive trends of |
---|
| 194 | !! temperature (idem for salinity) is given by: |
---|
| 195 | !! difft = 1/(e1t*e2t*e3t) { di-1[ aht e2u*e3u/e1u di(tb) ] |
---|
| 196 | !! + dj-1[ aht e1v*e3v/e2v dj(tb) ] } |
---|
| 197 | !! Note: key_zco defined, the e3t=e3u=e3v, the trend becomes: |
---|
| 198 | !! difft = 1/(e1t*e2t) { di-1[ aht e2u/e1u di(tb) ] |
---|
| 199 | !! + dj-1[ aht e1v/e2v dj(tb) ] } |
---|
| 200 | !! Add this trend to the general tracer trend (ta,sa): |
---|
| 201 | !! (ta,sa) = (ta,sa) + ( difft , diffs ) |
---|
| 202 | !! |
---|
| 203 | !! ** Action : - Update (ta,sa) arrays with the before iso-level |
---|
| 204 | !! harmonic mixing trend. |
---|
| 205 | !! |
---|
| 206 | !! History : |
---|
| 207 | !! 1.0 ! 87-06 (P. Andrich, D. L Hostis) Original code |
---|
| 208 | !! ! 91-11 (G. Madec) |
---|
| 209 | !! ! 95-11 (G. Madec) suppress volumetric scale factors |
---|
| 210 | !! ! 96-01 (G. Madec) statement function for e3 |
---|
| 211 | !! 8.5 ! 02-06 (G. Madec) F90: Free form and module |
---|
| 212 | !! 9.0 ! 04-08 (C. Talandier) New trends organization |
---|
| 213 | !! ! 05-11 (G. Madec) add zps case |
---|
| 214 | !! History of the tangent routine |
---|
| 215 | !! 9.0 ! 03-09 (F. Vigilant) adjoint of 9.0 |
---|
| 216 | !!---------------------------------------------------------------------- |
---|
| 217 | |
---|
| 218 | !! * Arguments |
---|
| 219 | INTEGER, INTENT( in ) :: kt ! ocean time-step index |
---|
| 220 | INTEGER , INTENT(in ) :: kit000 ! first time step index |
---|
| 221 | CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) |
---|
| 222 | INTEGER , INTENT(in ) :: kjpt ! number of tracers |
---|
| 223 | REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(inout) :: pgu_ad, pgv_ad ! tracer gradient at pstep levels |
---|
| 224 | REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb_ad ! before and now tracer fields |
---|
| 225 | REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta_ad ! tracer trend |
---|
| 226 | !! * Local declarations |
---|
| 227 | INTEGER :: ji, jj, jk, jn ! dummy loop indices |
---|
| 228 | INTEGER :: iku, ikv, ierr ! temporary integers |
---|
| 229 | REAL(wp) :: zabe1, zabe2, zbtr ! temporary scalars |
---|
| 230 | REAL(wp) :: ztaad, zsaad ! temporary scalars |
---|
| 231 | REAL(wp), POINTER, DIMENSION(:,:,:) :: ztuad, ztvad ! 3D workspace |
---|
| 232 | !!---------------------------------------------------------------------- |
---|
| 233 | ! |
---|
| 234 | CALL timing_start('tra_ldf_lap_adj') |
---|
| 235 | ! |
---|
| 236 | CALL wrk_alloc(jpi,jpj,jpk, ztuad, ztvad) |
---|
| 237 | ! |
---|
| 238 | ztvad(:,:,:) = 0.0_wp |
---|
| 239 | ztuad(:,:,:) = 0.0_wp |
---|
| 240 | ztaad = 0.0_wp |
---|
| 241 | zsaad = 0.0_wp |
---|
| 242 | |
---|
| 243 | IF( kt == nitend ) THEN |
---|
| 244 | IF(lwp) WRITE(numout,*) |
---|
| 245 | IF(lwp) WRITE(numout,*) 'tra_ldf_lap_adj : iso-level laplacian diffusion on ', cdtype |
---|
| 246 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' |
---|
| 247 | IF( .NOT. ALLOCATED( e1ur ) ) THEN |
---|
| 248 | ! This routine may be called for both active and passive tracers. |
---|
| 249 | ! Allocate and set saved arrays on first call only. |
---|
| 250 | ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr ) |
---|
| 251 | IF( lk_mpp ) CALL mpp_sum( ierr ) |
---|
| 252 | IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' ) |
---|
| 253 | ! |
---|
| 254 | e1ur(:,:) = e2u(:,:) / e1u(:,:) |
---|
| 255 | e2vr(:,:) = e1v(:,:) / e2v(:,:) |
---|
| 256 | ENDIF |
---|
| 257 | ENDIF |
---|
| 258 | |
---|
| 259 | |
---|
| 260 | DO jn = 1, kjpt |
---|
| 261 | ! ! ============= |
---|
| 262 | DO jk = 1, jpkm1 ! Vertical slab |
---|
| 263 | ! ! ============= |
---|
| 264 | |
---|
| 265 | ! 2. Second derivative (divergence) |
---|
| 266 | ! -------------------- |
---|
| 267 | DO jj = jpjm1, 2, -1 |
---|
| 268 | DO ji = fs_jpim1, fs_2, -1 ! vector opt. |
---|
| 269 | zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) ) |
---|
| 270 | ! horizontal diffusive trends |
---|
| 271 | ztaad = zbtr * pta_ad(ji,jj,jk,jn) |
---|
| 272 | ztuad(ji ,jj ,jk) = ztuad(ji ,jj ,jk) + ztaad |
---|
| 273 | ztuad(ji-1,jj ,jk) = ztuad(ji-1,jj ,jk) - ztaad |
---|
| 274 | ztvad(ji ,jj ,jk) = ztvad(ji ,jj ,jk) + ztaad |
---|
| 275 | ztvad(ji ,jj-1,jk) = ztvad(ji ,jj-1,jk) - ztaad |
---|
| 276 | !ztaad = 0.0_wp |
---|
| 277 | !zsaad = 0.0_wp |
---|
| 278 | END DO |
---|
| 279 | END DO |
---|
| 280 | |
---|
| 281 | ! 1. First derivative (gradient) |
---|
| 282 | ! ------------------- |
---|
| 283 | |
---|
| 284 | IF( ln_zps ) THEN ! set gradient at partial step level |
---|
| 285 | DO jj = jpjm1, 1, -1 |
---|
| 286 | DO ji = fs_jpim1, 1, -1 ! vector opt.? |
---|
| 287 | ! last level |
---|
| 288 | iku = mbku(ji,jj) |
---|
| 289 | ikv = mbkv(ji,jj) |
---|
| 290 | IF( iku == jk ) THEN |
---|
| 291 | zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku) |
---|
| 292 | pgu_ad(ji,jj,jn) = zabe1 * ztuad(ji,jj,jk) + pgu_ad(ji,jj,jn) |
---|
| 293 | ztuad(ji,jj,jk) = 0.0_wp |
---|
| 294 | ENDIF |
---|
| 295 | IF( ikv == jk ) THEN |
---|
| 296 | zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv) |
---|
| 297 | pgv_ad(ji,jj,jn) = zabe2 * ztvad(ji,jj,jk) + pgv_ad(ji,jj,jn) |
---|
| 298 | ztvad(ji,jj,jk) = 0.0_wp |
---|
| 299 | ENDIF |
---|
| 300 | END DO |
---|
| 301 | END DO |
---|
| 302 | ENDIF |
---|
| 303 | DO jj = jpjm1, 1, -1 |
---|
| 304 | DO ji = fs_jpim1, 1, -1 ! vector opt. ? |
---|
| 305 | zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk) |
---|
| 306 | zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk) |
---|
| 307 | ptb_ad(ji ,jj ,jk,jn) = ptb_ad(ji ,jj ,jk,jn) - (zabe1 * ztuad(ji,jj,jk) + zabe2 * ztvad(ji,jj,jk)) |
---|
| 308 | ptb_ad(ji+1,jj ,jk,jn) = ptb_ad(ji+1,jj ,jk,jn) + zabe1 * ztuad(ji,jj,jk) |
---|
| 309 | ptb_ad(ji ,jj+1,jk,jn) = ptb_ad(ji ,jj+1,jk,jn) + zabe2 * ztvad(ji,jj,jk) |
---|
| 310 | ztuad(ji,jj,jk) = 0.0_wp |
---|
| 311 | ztvad(ji,jj,jk) = 0.0_wp |
---|
| 312 | END DO |
---|
| 313 | END DO |
---|
| 314 | ! ! ============= |
---|
| 315 | END DO ! End of slab |
---|
| 316 | ! ! ============= |
---|
| 317 | END DO |
---|
| 318 | ! |
---|
| 319 | CALL timing_stop('tra_ldf_lap_adj') |
---|
| 320 | ! |
---|
| 321 | CALL wrk_dealloc(jpi,jpj,jpk, ztuad, ztvad) |
---|
| 322 | ! |
---|
| 323 | END SUBROUTINE tra_ldf_lap_adj |
---|
| 324 | |
---|
| 325 | |
---|
| 326 | SUBROUTINE tra_ldf_lap_adj_tst ( kumadt ) |
---|
| 327 | !!----------------------------------------------------------------------- |
---|
| 328 | !! |
---|
| 329 | !! *** ROUTINE example_adj_tst *** |
---|
| 330 | !! |
---|
| 331 | !! ** Purpose : Test the adjoint routine. |
---|
| 332 | !! |
---|
| 333 | !! ** Method : Verify the scalar product |
---|
| 334 | !! |
---|
| 335 | !! ( L dx )^T W dy = dx^T L^T W dy |
---|
| 336 | !! |
---|
| 337 | !! where L = tangent routine |
---|
| 338 | !! L^T = adjoint routine |
---|
| 339 | !! W = diagonal matrix of scale factors |
---|
| 340 | !! dx = input perturbation (random field) |
---|
| 341 | !! dy = L dx |
---|
| 342 | !! |
---|
| 343 | !! History : |
---|
| 344 | !! ! 08-08 (A. Vidard) |
---|
| 345 | !!----------------------------------------------------------------------- |
---|
| 346 | !! * Modules used |
---|
| 347 | |
---|
| 348 | !! * Arguments |
---|
| 349 | INTEGER, INTENT(IN) :: & |
---|
| 350 | & kumadt ! Output unit |
---|
| 351 | |
---|
| 352 | !! * Local declarations |
---|
| 353 | INTEGER :: & |
---|
| 354 | & ji, & ! dummy loop indices |
---|
| 355 | & jj, & |
---|
| 356 | & jk |
---|
| 357 | INTEGER, DIMENSION(jpi,jpj) :: & |
---|
| 358 | & iseed_2d ! 2D seed for the random number generator |
---|
| 359 | REAL(KIND=wp) :: & |
---|
| 360 | & zsp1, & ! scalar product involving the tangent routine |
---|
| 361 | & zsp1_T, & |
---|
| 362 | & zsp1_S, & |
---|
| 363 | & zsp2, & ! scalar product involving the adjoint routine |
---|
| 364 | & zsp2_1, & |
---|
| 365 | & zsp2_2, & |
---|
| 366 | & zsp2_3, & |
---|
| 367 | & zsp2_4, & |
---|
| 368 | & zsp2_5, & |
---|
| 369 | & zsp2_6, & |
---|
| 370 | & zsp2_7, & |
---|
| 371 | & zsp2_8, & |
---|
| 372 | & zsp2_T, & |
---|
| 373 | & zsp2_S |
---|
| 374 | REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
| 375 | & ztb_tlin , & ! Tangent input |
---|
| 376 | & zsb_tlin , & ! Tangent input |
---|
| 377 | & zta_tlin , & ! Tangent input |
---|
| 378 | & zsa_tlin , & ! Tangent input |
---|
| 379 | & zta_tlout, & ! Tangent output |
---|
| 380 | & zsa_tlout, & ! Tangent output |
---|
| 381 | & zta_adin, & ! Adjoint input |
---|
| 382 | & zsa_adin, & ! Adjoint input |
---|
| 383 | & ztb_adout , & ! Adjoint output |
---|
| 384 | & zsb_adout , & ! Adjoint output |
---|
| 385 | & zta_adout , & ! Adjoint output |
---|
| 386 | & zsa_adout , & ! Adjoint output |
---|
| 387 | & z3r ! 3D random field |
---|
| 388 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 389 | & zgtu_tlin , & ! Tangent input |
---|
| 390 | & zgsu_tlin , & ! Tangent input |
---|
| 391 | & zgtv_tlin , & ! Tangent input |
---|
| 392 | & zgsv_tlin , & ! Tangent input |
---|
| 393 | & zgtu_adout , & ! Adjoint output |
---|
| 394 | & zgsu_adout , & ! Adjoint output |
---|
| 395 | & zgtv_adout , & ! Adjoint output |
---|
| 396 | & zgsv_adout , & ! Adjoint output |
---|
| 397 | & z2r ! 2D random field |
---|
| 398 | CHARACTER(LEN=14) :: cl_name |
---|
| 399 | ! Allocate memory |
---|
| 400 | |
---|
| 401 | ALLOCATE( & |
---|
| 402 | & ztb_tlin(jpi,jpj,jpk), & |
---|
| 403 | & zsb_tlin(jpi,jpj,jpk), & |
---|
| 404 | & zta_tlin(jpi,jpj,jpk), & |
---|
| 405 | & zsa_tlin(jpi,jpj,jpk), & |
---|
| 406 | & zgtu_tlin(jpi,jpj), & |
---|
| 407 | & zgsu_tlin(jpi,jpj), & |
---|
| 408 | & zgtv_tlin(jpi,jpj), & |
---|
| 409 | & zgsv_tlin(jpi,jpj), & |
---|
| 410 | & zta_tlout(jpi,jpj,jpk), & |
---|
| 411 | & zsa_tlout(jpi,jpj,jpk), & |
---|
| 412 | & zta_adin(jpi,jpj,jpk), & |
---|
| 413 | & zsa_adin(jpi,jpj,jpk), & |
---|
| 414 | & ztb_adout(jpi,jpj,jpk), & |
---|
| 415 | & zsb_adout(jpi,jpj,jpk), & |
---|
| 416 | & zta_adout(jpi,jpj,jpk), & |
---|
| 417 | & zsa_adout(jpi,jpj,jpk), & |
---|
| 418 | & zgtu_adout(jpi,jpj), & |
---|
| 419 | & zgsu_adout(jpi,jpj), & |
---|
| 420 | & zgtv_adout(jpi,jpj), & |
---|
| 421 | & zgsv_adout(jpi,jpj), & |
---|
| 422 | & z3r(jpi,jpj,jpk), & |
---|
| 423 | & z2r(jpi,jpj) & |
---|
| 424 | & ) |
---|
| 425 | |
---|
| 426 | |
---|
| 427 | !======================================================================= |
---|
| 428 | ! 1) dx = ( tb_tl, ta_tl, sb_tl, sa_tl, gtu_tl, gtv_tl, gsu_tl, gsv_tl ) |
---|
| 429 | ! dy = ( ta_tl, sa_tl ) |
---|
| 430 | !======================================================================= |
---|
| 431 | |
---|
| 432 | !-------------------------------------------------------------------- |
---|
| 433 | ! Reset the tangent and adjoint variables |
---|
| 434 | !-------------------------------------------------------------------- |
---|
| 435 | |
---|
| 436 | ztb_tlin(:,:,:) = 0.0_wp |
---|
| 437 | zsb_tlin(:,:,:) = 0.0_wp |
---|
| 438 | zta_tlin(:,:,:) = 0.0_wp |
---|
| 439 | zsa_tlin(:,:,:) = 0.0_wp |
---|
| 440 | zgtu_tlin(:,:) = 0.0_wp |
---|
| 441 | zgsu_tlin(:,:) = 0.0_wp |
---|
| 442 | zgtv_tlin(:,:) = 0.0_wp |
---|
| 443 | zgsv_tlin(:,:) = 0.0_wp |
---|
| 444 | zta_tlout(:,:,:) = 0.0_wp |
---|
| 445 | zsa_tlout(:,:,:) = 0.0_wp |
---|
| 446 | zta_adin(:,:,:) = 0.0_wp |
---|
| 447 | zsa_adin(:,:,:) = 0.0_wp |
---|
| 448 | ztb_adout(:,:,:) = 0.0_wp |
---|
| 449 | zsb_adout(:,:,:) = 0.0_wp |
---|
| 450 | zta_adout(:,:,:) = 0.0_wp |
---|
| 451 | zsa_adout(:,:,:) = 0.0_wp |
---|
| 452 | zgtu_adout(:,:) = 0.0_wp |
---|
| 453 | zgsu_adout(:,:) = 0.0_wp |
---|
| 454 | zgtv_adout(:,:) = 0.0_wp |
---|
| 455 | zgsv_adout(:,:) = 0.0_wp |
---|
| 456 | |
---|
| 457 | tsb_tl(:,:,:,:) = 0.0_wp |
---|
| 458 | tsa_tl(:,:,:,:) = 0.0_wp |
---|
| 459 | gtsu_tl(:,:,:) = 0.0_wp |
---|
| 460 | gtsv_tl(:,:,:) = 0.0_wp |
---|
| 461 | tsb_ad(:,:,:,:) = 0.0_wp |
---|
| 462 | tsa_ad(:,:,:,:) = 0.0_wp |
---|
| 463 | gtsu_ad(:,:,:) = 0.0_wp |
---|
| 464 | gtsv_ad(:,:,:) = 0.0_wp |
---|
| 465 | |
---|
| 466 | !-------------------------------------------------------------------- |
---|
| 467 | ! Initialize the tangent input with random noise: dx |
---|
| 468 | !-------------------------------------------------------------------- |
---|
| 469 | CALL grid_random( z3r, 'T', 0.0_wp, stdt ) |
---|
| 470 | ztb_tlin(:,:,:) = z3r(:,:,:) |
---|
| 471 | CALL grid_random( z3r, 'T', 0.0_wp, stds ) |
---|
| 472 | zsb_tlin(:,:,:) = z3r(:,:,:) |
---|
| 473 | CALL grid_random( z3r, 'T', 0.0_wp, stdt ) |
---|
| 474 | zta_tlin(:,:,:) = z3r(:,:,:) |
---|
| 475 | CALL grid_random( z3r, 'T', 0.0_wp, stds ) |
---|
| 476 | zsa_tlin(:,:,:) = z3r(:,:,:) |
---|
| 477 | CALL grid_random( z2r, 'U', 0.0_wp, stds ) |
---|
| 478 | zgtu_tlin(:,:) = z2r(:,:) |
---|
| 479 | CALL grid_random( z2r, 'U', 0.0_wp, stds ) |
---|
| 480 | zgsu_tlin(:,:) = z2r(:,:) |
---|
| 481 | CALL grid_random( z2r, 'V', 0.0_wp, stds ) |
---|
| 482 | zgtv_tlin(:,:) = z2r(:,:) |
---|
| 483 | CALL grid_random( z2r, 'V', 0.0_wp, stds ) |
---|
| 484 | zgsv_tlin(:,:) = z2r(:,:) |
---|
| 485 | |
---|
| 486 | tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:) |
---|
| 487 | tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:) |
---|
| 488 | tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:) |
---|
| 489 | tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:) |
---|
| 490 | gtsu_tl(:,:,jp_tem) = zgtu_tlin(:,:) |
---|
| 491 | gtsu_tl(:,:,jp_sal) = zgsu_tlin(:,:) |
---|
| 492 | gtsv_tl(:,:,jp_tem) = zgtv_tlin(:,:) |
---|
| 493 | gtsv_tl(:,:,jp_sal) = zgsv_tlin(:,:) |
---|
| 494 | |
---|
| 495 | CALL tra_ldf_lap_tan( nit000, nit000, 'TRA', gtsu_tl, gtsv_tl, tsb_tl, tsa_tl, jpts ) |
---|
| 496 | |
---|
| 497 | zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem) |
---|
| 498 | zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal) |
---|
| 499 | |
---|
| 500 | !-------------------------------------------------------------------- |
---|
| 501 | ! Initialize the adjoint variables: dy^* = W dy |
---|
| 502 | !-------------------------------------------------------------------- |
---|
| 503 | |
---|
| 504 | DO jk = 1, jpk |
---|
| 505 | DO jj = nldj, nlej |
---|
| 506 | DO ji = nldi, nlei |
---|
| 507 | zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) & |
---|
| 508 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
| 509 | & * tmask(ji,jj,jk) * wesp_s(jk) |
---|
| 510 | zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) & |
---|
| 511 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
| 512 | & * tmask(ji,jj,jk) * wesp_t(jk) |
---|
| 513 | END DO |
---|
| 514 | END DO |
---|
| 515 | END DO |
---|
| 516 | |
---|
| 517 | !-------------------------------------------------------------------- |
---|
| 518 | ! Compute the scalar product: ( L dx )^T W dy |
---|
| 519 | !------------------------------------------------------------------- |
---|
| 520 | |
---|
| 521 | zsp1_T = DOT_PRODUCT( zta_tlout, zta_adin ) |
---|
| 522 | zsp1_S = DOT_PRODUCT( zsa_tlout, zsa_adin ) |
---|
| 523 | zsp1 = zsp1_T + zsp1_S |
---|
| 524 | |
---|
| 525 | !-------------------------------------------------------------------- |
---|
| 526 | ! Call the adjoint routine: dx^* = L^T dy^* |
---|
| 527 | !-------------------------------------------------------------------- |
---|
| 528 | |
---|
| 529 | tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:) |
---|
| 530 | tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:) |
---|
| 531 | |
---|
| 532 | CALL tra_ldf_lap_adj( nit000, nit000, 'TRA', gtsu_ad, gtsv_ad, tsb_ad, tsa_ad, jpts ) |
---|
| 533 | |
---|
| 534 | ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem) |
---|
| 535 | zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal) |
---|
| 536 | zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem) |
---|
| 537 | zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal) |
---|
| 538 | zgtu_adout(:,:) = gtsu_ad(:,:,jp_tem) |
---|
| 539 | zgsu_adout(:,:) = gtsu_ad(:,:,jp_sal) |
---|
| 540 | zgtv_adout(:,:) = gtsv_ad(:,:,jp_tem) |
---|
| 541 | zgsv_adout(:,:) = gtsv_ad(:,:,jp_sal) |
---|
| 542 | |
---|
| 543 | zsp2_1 = DOT_PRODUCT( ztb_tlin , ztb_adout ) |
---|
| 544 | zsp2_2 = DOT_PRODUCT( zta_tlin , zta_adout ) |
---|
| 545 | zsp2_3 = DOT_PRODUCT( zgtu_tlin, zgtu_adout ) |
---|
| 546 | zsp2_4 = DOT_PRODUCT( zgtv_tlin, zgtv_adout ) |
---|
| 547 | zsp2_5 = DOT_PRODUCT( zsb_tlin , zsb_adout ) |
---|
| 548 | zsp2_6 = DOT_PRODUCT( zsa_tlin , zsa_adout ) |
---|
| 549 | zsp2_7 = DOT_PRODUCT( zgsu_tlin, zgsu_adout ) |
---|
| 550 | zsp2_8 = DOT_PRODUCT( zgsv_tlin, zgsv_adout ) |
---|
| 551 | |
---|
| 552 | zsp2_T = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 |
---|
| 553 | zsp2_S = zsp2_5 + zsp2_6 + zsp2_7 + zsp2_8 |
---|
| 554 | zsp2 = zsp2_T + zsp2_S |
---|
| 555 | |
---|
| 556 | cl_name = 'tra_ldf_lap_ad' |
---|
| 557 | CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) |
---|
| 558 | |
---|
| 559 | DEALLOCATE( & |
---|
| 560 | & ztb_tlin, & ! Tangent input |
---|
| 561 | & zsb_tlin, & ! Tangent input |
---|
| 562 | & zta_tlin, & ! Tangent input |
---|
| 563 | & zsa_tlin, & ! Tangent input |
---|
| 564 | & zgtu_tlin, & ! Tangent input |
---|
| 565 | & zgsu_tlin, & ! Tangent input |
---|
| 566 | & zgtv_tlin, & ! Tangent input |
---|
| 567 | & zgsv_tlin, & ! Tangent input |
---|
| 568 | & zta_tlout, & ! Tangent output |
---|
| 569 | & zsa_tlout, & ! Tangent output |
---|
| 570 | & zta_adin, & ! Adjoint input |
---|
| 571 | & zsa_adin, & ! Adjoint input |
---|
| 572 | & ztb_adout, & ! Adjoint output |
---|
| 573 | & zsb_adout, & ! Adjoint output |
---|
| 574 | & zta_adout, & ! Adjoint output |
---|
| 575 | & zsa_adout, & ! Adjoint output |
---|
| 576 | & zgtu_adout, & ! Adjoint output |
---|
| 577 | & zgsu_adout, & ! Adjoint output |
---|
| 578 | & zgtv_adout, & ! Adjoint output |
---|
| 579 | & zgsv_adout, & ! Adjoint output |
---|
| 580 | & z3r, & ! 3D random field |
---|
| 581 | & z2r & |
---|
| 582 | & ) |
---|
| 583 | |
---|
| 584 | |
---|
| 585 | END SUBROUTINE tra_ldf_lap_adj_tst |
---|
| 586 | #endif |
---|
| 587 | !!============================================================================== |
---|
| 588 | END MODULE traldf_lap_tam |
---|