Changeset 455 for trunk/NEMO
- Timestamp:
- 2006-05-10T18:53:54+02:00 (18 years ago)
- Location:
- trunk/NEMO/OPA_SRC/DYN
- Files:
-
- 3 deleted
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/divcur.F90
r392 r455 51 51 !! (note that the Asselin filter has not been applied on hdivb) 52 52 !! - compute the now divergence given by : 53 !! * s-coordinate ('key_s_coord' defined)54 53 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 55 !! * z-coordinate (default key)56 !! hdivn = 1/(e1t*e2t) [ di(e2u un) + dj(e1v vn) ]54 !! Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the 55 !! above expression 57 56 !! - apply lateral boundary conditions on hdivn 58 57 !! II. vorticity : … … 109 108 DO jj = 2, jpjm1 110 109 DO ji = fs_2, fs_jpim1 ! vector opt. 111 #if defined key_s_coord || defined key_partial_steps 112 hdivn(ji,jj,jk) = & 113 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 114 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 115 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 116 #else 110 #if defined key_zco 117 111 hdivn(ji,jj,jk) = ( e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj ) * un(ji-1,jj ,jk) & 118 112 & + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji ,jj-1) * vn(ji ,jj-1,jk) ) & 119 & / ( e1t(ji,jj) * e2t(ji,jj) ) 113 & / ( e1t(ji,jj) * e2t(ji,jj) ) 114 #else 115 hdivn(ji,jj,jk) = & 116 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 117 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 118 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 120 119 #endif 121 120 END DO … … 130 129 IF( lp_obc_south ) hdivn(nis0 :nis1 ,njs0 :njs1 ,jk) = 0.e0 ! south 131 130 #endif 132 #if defined key_agrif133 if ( .NOT. AGRIF_Root() ) then134 IF ((nbondi == 1).OR.(nbondi == 2)) hdivn(nlci-1 , : ,jk) = 0.e0 ! east135 IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2 , : ,jk) = 0.e0 ! west136 IF ((nbondj == 1).OR.(nbondj == 2)) hdivn(: ,nlcj-1 ,jk) = 0.e0 ! north137 IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(: ,2 ,jk) = 0.e0 ! south138 endif139 #endif140 131 141 132 ! ! -------- … … 260 251 !! (note that the Asselin filter has not been applied on hdivb) 261 252 !! - compute the now divergence given by : 262 !! * s-coordinate ('key_s_coord' defined)263 253 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 264 !! * z-coordinate (default key)265 !! hdivn = 1/(e1t*e2t) [ di(e2u un) + dj(e1v vn) ]254 !! Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the 255 !! above expression 266 256 !! - apply lateral boundary conditions on hdivn 267 257 !! - Relavtive Vorticity : … … 313 303 DO jj = 2, jpjm1 314 304 DO ji = fs_2, fs_jpim1 ! vector opt. 315 #if defined key_s_coord || defined key_partial_steps 316 hdivn(ji,jj,jk) = & 317 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 318 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 319 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 320 #else 305 #if defined key_zco 321 306 hdivn(ji,jj,jk) = ( e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj ) * un(ji-1,jj ,jk) & 322 307 & + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji ,jj-1) * vn(ji ,jj-1,jk) ) & 323 308 / ( e1t(ji,jj) * e2t(ji,jj) ) 309 #else 310 hdivn(ji,jj,jk) = & 311 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 312 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 313 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 324 314 #endif 325 315 END DO … … 334 324 IF( lp_obc_south ) hdivn(nis0 :nis1 ,njs0 :njs1 ,jk) = 0.e0 ! south 335 325 #endif 336 #if defined key_agrif337 if ( .NOT. AGRIF_Root() ) then338 IF ((nbondi == 1).OR.(nbondi == 2)) hdivn(nlci-1 , : ,jk) = 0.e0 ! east339 IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2 , : ,jk) = 0.e0 ! west340 IF ((nbondj == 1).OR.(nbondj == 2)) hdivn(: ,nlcj-1 ,jk) = 0.e0 ! north341 IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(: ,2 ,jk) = 0.e0 ! south342 endif343 #endif344 326 ! ! -------- 345 327 ! relative vorticity ! rot -
trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
r359 r455 6 6 7 7 !!---------------------------------------------------------------------- 8 !! dyn_hpg : update the momentum trend with the horizontal8 !! dyn_hpg : update the momentum trend with the now horizontal 9 9 !! gradient of the hydrostatic pressure 10 !! 11 !! default case : use of 3D work arrays (vector opt. available) 12 !! key_s_coord : s-coordinate 13 !! key_partial_steps : z-coordinate with partial steps 14 !! default key : z-coordinate 10 !! default case : k-j-i loops (vector opt. available) 11 !! hpg_ctl : initialisation and control of options 12 !! hpg_zco : z-coordinate scheme 13 !! hpg_zps : z-coordinate plus partial steps (interpolation) 14 !! hpg_sco : s-coordinate (standard jacobian formulation) 15 !! hpg_hel : s-coordinate (helsinki modification) 16 !! hpg_wdj : s-coordinate (weighted density jacobian) 17 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 18 !! hpg_rot : s-coordinate (ROTated axes scheme) 15 19 !!---------------------------------------------------------------------- 16 20 !! * Modules used 17 21 USE oce ! ocean dynamics and tracers 18 22 USE dom_oce ! ocean space and time domain 23 USE dynhpg_jki ! 19 24 USE phycst ! physical constants 20 25 USE in_out_manager ! I/O manager … … 22 27 USE trdmod_oce ! ocean variables trends 23 28 USE prtctl ! Print control 29 USE lbclnk ! lateral boundary condition 24 30 25 31 IMPLICIT NONE … … 28 34 !! * Accessibility 29 35 PUBLIC dyn_hpg ! routine called by step.F90 36 37 #if defined key_mpp_omp 38 !!---------------------------------------------------------------------- 39 !! 'key_mpp_omp' : j-k-i loop (j-slab) 40 !!---------------------------------------------------------------------- 41 LOGICAL, PUBLIC, PARAMETER :: lk_dynhpg_jki = .TRUE. !: OpenMP hpg flag 42 LOGICAL, PUBLIC, PARAMETER :: lk_dynhpg = .FALSE. !: vector hpg flag 43 #else 44 !!---------------------------------------------------------------------- 45 !! default case : k-j-i loop (vector opt.) 46 !!---------------------------------------------------------------------- 47 LOGICAL, PUBLIC, PARAMETER :: lk_dynhpg_jki = .FALSE. !: OpenMP hpg flag 48 LOGICAL, PUBLIC, PARAMETER :: lk_dynhpg = .TRUE. !: vector hpg flag 49 #endif 50 51 !! * Share module variables 52 LOGICAL :: & !!! ** nam_dynhpg ** hpg flags 53 ln_hpg_zco = .TRUE. , & ! z-coordinate - full steps 54 ln_hpg_zps = .FALSE., & ! z-coordinate - partial steps (interpolation) 55 ln_hpg_sco = .FALSE., & ! s-coordinate (standard jacobian formulation) 56 ln_hpg_hel = .FALSE., & ! s-coordinate (helsinki modification) 57 ln_hpg_wdj = .FALSE., & ! s-coordinate (weighted density jacobian) 58 ln_hpg_djc = .FALSE., & ! s-coordinate (Density Jacobian with Cubic polynomial) 59 ln_hpg_rot = .FALSE. ! s-coordinate (ROTated axes scheme) 60 61 REAL(wp) :: & !!! ** nam_dynhpg ** 62 gamm = 0.e0 ! weighting coefficient 63 64 INTEGER :: & ! 65 nhpg = 0 ! = 0 to 6, type of pressure gradient scheme used 66 ! ! (deduced from ln_hpg_... flags) 30 67 31 68 !! * Substitutions … … 40 77 CONTAINS 41 78 42 #if defined key_s_coord43 !!----------------------------------------------------------------------44 !! 'key_s_coord' : s-coordinate45 !!----------------------------------------------------------------------46 47 79 SUBROUTINE dyn_hpg( kt ) 48 80 !!--------------------------------------------------------------------- 49 81 !! *** ROUTINE dyn_hpg *** 50 82 !! 51 !! ** Purpose : Compute the now momentum trend due to the hor. gradient 52 !! of the hydrostatic pressure. Add it to the general momentum trend. 53 !! 54 !! ** Method : The now hydrostatic pressure gradient at a given level 55 !! jk is computed by taking the vertical integral of the in-situ 56 !! density gradient along the model level from the suface to that 57 !! level. s-coordinates ('key_s_coord'): a corrective term is added 58 !! to the horizontal pressure gradient : 59 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 60 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 83 !! ** Method : Call the hydrostatic pressure gradient routine 84 !! using the scheme defined in the namelist (nhpg parameter) 85 !! 86 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 87 !! - Save the trend (l_trddyn=T) 88 !! - Control print (ln_ctl) 89 !! 90 !! History : 91 !! 9.0 ! 05-10 (A. Beckmann, G. Madec) various s-coordinate options 92 !!---------------------------------------------------------------------- 93 !! * Arguments 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 95 96 !! * local declarations 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 98 ztrdu, ztrdv ! 3D temporary workspace 99 !!---------------------------------------------------------------------- 100 101 IF( kt == nit000 ) CALL hpg_ctl ! initialisation & control of options 102 103 ! Temporary saving of ua and va trends (l_trddyn) 104 IF( l_trddyn ) THEN 105 ztrdu(:,:,:) = ua(:,:,:) 106 ztrdv(:,:,:) = va(:,:,:) 107 ENDIF 108 109 SELECT CASE ( nhpg ) ! Hydrastatic pressure gradient computation 110 CASE ( 0 ) ! z-coordinate 111 CALL hpg_zco( kt ) 112 CASE ( 1 ) ! z-coordinate plus partial steps (interpolation) 113 CALL hpg_zps( kt ) 114 CASE ( 2 ) ! s-coordinate (standard jacobian formulation) 115 CALL hpg_sco( kt ) 116 CASE ( 3 ) ! s-coordinate (helsinki modification) 117 CALL hpg_hel( kt ) 118 CASE ( 4 ) ! s-coordinate (weighted density jacobian) 119 CALL hpg_wdj( kt ) 120 CASE ( 5 ) ! s-coordinate (Density Jacobian with Cubic polynomial) 121 CALL hpg_djc( kt ) 122 CASE ( 6 ) ! s-coordinate (ROTated axes scheme) 123 CALL hpg_rot( kt ) 124 CASE ( 10 ) ! z-coordinate 125 CALL hpg_zco_jki( kt ) 126 CASE ( 11 ) ! z-coordinate plus partial steps (interpolation) 127 CALL hpg_zps_jki( kt ) 128 CASE ( 12 ) ! s-coordinate (standard jacobian formulation) 129 CALL hpg_sco_jki( kt ) 130 END SELECT 131 132 ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 133 IF( l_trddyn ) THEN 134 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 135 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 136 CALL trd_mod( ztrdu, ztrdv, jpdtdhpg, 'DYN', kt ) 137 ENDIF 138 139 IF(ln_ctl) THEN ! print sum trends (used for debugging) 140 CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, & 141 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 142 ENDIF 143 144 END SUBROUTINE dyn_hpg 145 146 147 SUBROUTINE hpg_ctl 148 !!---------------------------------------------------------------------- 149 !! *** ROUTINE hpg_ctl *** 150 !! 151 !! ** Purpose : initializations for the hydrostatic pressure gradient 152 !! computation and consistency control 153 !! 154 !! ** Action : Read the namelist namdynhpg and check the consistency 155 !! with the type of vertical coordinate used (zco, zps, sco) 156 !! 157 !! History : 158 !! 9.0 ! 05-10 (A. Beckmann) Original code 159 !!---------------------------------------------------------------------- 160 INTEGER :: ioptio = 0 ! temporary integer 161 162 NAMELIST/nam_dynhpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 163 & ln_hpg_hel, ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, & 164 & gamm 165 !!---------------------------------------------------------------------- 166 167 ! Read Namelist nam_dynhpg : pressure gradient calculation options 168 REWIND ( numnam ) 169 READ ( numnam, nam_dynhpg ) 170 171 ! Control print 172 IF(lwp) THEN 173 WRITE(numout,*) 174 WRITE(numout,*) 'dyn:hpg_ctl : hydrostatic pressure gradient control' 175 WRITE(numout,*) '~~~~~~~~~~~' 176 WRITE(numout,*) ' Namelist nam_dynhpg : choice of hpg scheme' 177 WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco 178 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 179 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 180 WRITE(numout,*) ' s-coord. (helsinki modification) ln_hpg_hel = ', ln_hpg_hel 181 WRITE(numout,*) ' s-coord. (weighted density jacobian) ln_hpg_wdj = ', ln_hpg_wdj 182 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 183 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 184 WRITE(numout,*) ' weighting coeff. (wdj scheme) gamm = ', gamm 185 ENDIF 186 187 ! set nhpg from ln_hpg_... flags 188 IF( ln_hpg_zco ) nhpg = 0 189 IF( ln_hpg_zps ) nhpg = 1 190 IF( ln_hpg_sco ) nhpg = 2 191 IF( ln_hpg_hel ) nhpg = 3 192 IF( ln_hpg_wdj ) nhpg = 4 193 IF( ln_hpg_djc ) nhpg = 5 194 IF( ln_hpg_rot ) nhpg = 6 195 196 ! Consitency check 197 ioptio = 0 198 IF( ln_hpg_zco ) ioptio = ioptio + 1 199 IF( ln_hpg_zps ) ioptio = ioptio + 1 200 IF( ln_hpg_sco ) ioptio = ioptio + 1 201 IF( ln_hpg_hel ) ioptio = ioptio + 1 202 IF( ln_hpg_wdj ) ioptio = ioptio + 1 203 IF( ln_hpg_djc ) ioptio = ioptio + 1 204 IF( ln_hpg_rot ) ioptio = ioptio + 1 205 IF ( ioptio > 1 ) THEN 206 IF(lwp) WRITE(numout,cform_err) 207 IF(lwp) WRITE(numout,*) ' several hydrostatic pressure gradient options used' 208 nstop = nstop + 1 209 ENDIF 210 211 IF( lk_dynhpg_jki ) THEN 212 nhpg = nhpg + 10 213 IF(lwp) WRITE(numout,*) 214 IF(lwp) WRITE(numout,*) ' Autotasking or OPENMP: use j-k-i loops (i.e. _jki routines)' 215 ENDIF 216 217 END SUBROUTINE hpg_ctl 218 219 220 SUBROUTINE hpg_zco( kt ) 221 !!--------------------------------------------------------------------- 222 !! *** ROUTINE hpg_zco *** 223 !! 224 !! ** Method : z-coordinate case, levels are horizontal surfaces. 225 !! The now hydrostatic pressure gradient at a given level, jk, 226 !! is computed by taking the vertical integral of the in-situ 227 !! density gradient along the model level from the suface to that 228 !! level: zhpi = grav ..... 229 !! zhpj = grav ..... 61 230 !! add it to the general momentum trend (ua,va). 62 !! ua = ua - 1/e1u * zhpi63 !! va = va - 1/e2v * zhpj64 !! 231 !! ua = ua - 1/e1u * zhpi 232 !! va = va - 1/e2v * zhpj 233 !! 65 234 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 66 !! - Save the trend in (utrd,vtrd) ('key_trddyn')67 235 !! 68 236 !! History : 69 !! 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code 70 !! ! 91-11 (G. Madec) 71 !! ! 96-01 (G. Madec) s-coordinates 72 !! ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg 73 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module, vector opt. 74 !! 9.0 ! 04-08 (C. Talandier) New trends organization 237 !! 1.0 ! 87-09 (P. Andrich, M.-A. Foujols) Original code 238 !! 5.0 ! 91-11 (G. Madec) 239 !! 7.0 ! 96-01 (G. Madec) 240 !! 8.0 ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg 241 !! 8.5 ! 02-07 (G. Madec) F90: Free form and module 75 242 !!---------------------------------------------------------------------- 76 243 !! * modules used 77 244 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 78 245 & zhpj => sa ! use sa as 3D workspace 79 246 80 247 !! * Arguments 81 248 INTEGER, INTENT( in ) :: kt ! ocean time-step index 82 249 83 !! * Local declarations250 !! * local declarations 84 251 INTEGER :: ji, jj, jk ! dummy loop indices 85 REAL(wp) :: & 86 zcoef0, zcoef1, zuap, zvap ! temporary scalars 87 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 88 ztdua, ztdva ! temporary scalars 89 !!---------------------------------------------------------------------- 90 252 REAL(wp) :: & 253 zcoef0, zcoef1 ! temporary scalars 254 !!---------------------------------------------------------------------- 255 91 256 IF( kt == nit000 ) THEN 92 257 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) 'dyn _hpg: hydrostatic pressure gradient trend'94 IF(lwp) WRITE(numout,*) '~~~~~~~ s-coordinate case, vector opt. case'258 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 259 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 95 260 ENDIF 96 97 ! Save ua and va trends 98 IF( l_trddyn ) THEN 99 ztdua(:,:,:) = ua(:,:,:) 100 ztdva(:,:,:) = va(:,:,:) 101 ENDIF 102 103 ! 0. Local constant initialization 104 ! -------------------------------- 261 262 263 ! Local constant initialization 264 ! ----------------------------- 105 265 zcoef0 = - grav * 0.5 106 zuap = 0.e0 107 zvap = 0.e0 108 109 ! 1. Surface value 110 ! ---------------- 111 DO jj = 2, jpjm1 112 DO ji = fs_2, fs_jpim1 ! vector opt. 113 ! hydrostatic pressure gradient along s-surfaces 114 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) & 115 * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1) ) 116 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) & 117 * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1) ) 118 ! s-coordinate pressure gradient correction 119 zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) ) & 120 * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 121 zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) ) & 122 * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 123 ! add to the general momentum trend 124 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 125 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 126 END DO 127 END DO 128 129 ! 2. interior value (2=<jk=<jpkm1) 130 ! ----------------- 131 DO jk = 2, jpkm1 132 DO jj = 2, jpjm1 133 DO ji = fs_2, fs_jpim1 ! vector opt. 134 ! hydrostatic pressure gradient along s-surfaces 135 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 136 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 137 & -fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) 138 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 139 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 140 & -fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) 141 ! s-coordinate pressure gradient correction 142 zuap = -zcoef0 * ( rhd(ji+1,jj ,jk) + rhd(ji,jj,jk) ) & 143 * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 144 zvap = -zcoef0 * ( rhd(ji ,jj+1,jk) + rhd(ji,jj,jk) ) & 145 * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 146 ! add to the general momentum trend 147 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 148 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 149 END DO 150 END DO 151 END DO 152 153 ! save the hydrostatic pressure gradient trends for diagnostic 154 ! momentum trends 155 IF( l_trddyn ) THEN 156 zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 157 zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 158 CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 159 ENDIF 160 161 IF(ln_ctl) THEN ! print sum trends (used for debugging) 162 CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, & 163 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 164 ENDIF 165 166 END SUBROUTINE dyn_hpg 167 168 #elif defined key_partial_steps 169 !!--------------------------------------------------------------------- 170 !! 'key_partial_steps' z-coordinate partial steps 171 !!--------------------------------------------------------------------- 172 173 SUBROUTINE dyn_hpg( kt ) 174 !!--------------------------------------------------------------------- 175 !! *** ROUTINE dyn_hpg *** 176 !! 177 !! ** Purpose : Compute the now momentum trend due to the horizontal 178 !! gradient of the hydrostatic pressure. Add it to the general 179 !! momentum trend. 180 !! 181 !! ** Method : The now hydrostatic pressure gradient at a given level 182 !! jk is computed by taking the vertical integral of the in-situ 183 !! density gradient along the model level from the suface to that 184 !! level: zhpi = grav ..... 185 !! zhpj = grav ..... 186 !! add it to the general momentum trend (ua,va). 187 !! ua = ua - 1/e1u * zhpi 188 !! va = va - 1/e2v * zhpj 189 !! 190 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 191 !! - Save the trend in (utrd,vtrd) ('key_trddyn') 192 !! 193 !! History : 194 !! 8.5 ! 02-08 (A. Bozec) Original code 195 !!---------------------------------------------------------------------- 196 !! * modules used 197 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 198 & zhpj => sa ! use sa as 3D workspace 199 200 !! * Arguments 201 INTEGER, INTENT( in ) :: kt ! ocean time-step index 202 203 !! * local declarations 204 INTEGER :: ji, jj, jk ! dummy loop indices 205 INTEGER :: iku, ikv ! temporary integers 206 REAL(wp) :: & 207 zcoef0, zcoef1, zuap, & ! temporary scalars 208 zcoef2, zcoef3, zvap ! " " 209 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 210 ztdua, ztdva ! temporary scalars 211 !!---------------------------------------------------------------------- 212 213 IF( kt == nit000 ) THEN 214 IF(lwp) WRITE(numout,*) 215 IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' 216 IF(lwp) WRITE(numout,*) '~~~~~~~ z-coordinate with partial steps' 217 IF(lwp) WRITE(numout,*) ' vector optimization, no autotasking' 218 ENDIF 219 220 ! Save ua and va trends 221 IF( l_trddyn ) THEN 222 ztdua(:,:,:) = ua(:,:,:) 223 ztdva(:,:,:) = va(:,:,:) 224 ENDIF 225 226 ! 0. Local constant initialization 227 ! -------------------------------- 228 zcoef0 = - grav * 0.5 229 zuap = 0.e0 230 zvap = 0.e0 231 232 ! 1. Surface value 233 ! ---------------- 266 267 ! Surface value 268 ! ------------- 234 269 DO jj = 2, jpjm1 235 270 DO ji = fs_2, fs_jpim1 ! vector opt. … … 244 279 END DO 245 280 281 ! interior value (2=<jk=<jpkm1) 282 ! -------------- 283 DO jk = 2, jpkm1 284 DO jj = 2, jpjm1 285 DO ji = fs_2, fs_jpim1 ! vector opt. 286 zcoef1 = zcoef0 * fse3w(ji,jj,jk) 287 ! hydrostatic pressure gradient 288 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 289 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 290 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj) 291 292 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 293 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 294 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj) 295 ! add to the general momentum trend 296 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 297 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 298 END DO 299 END DO 300 END DO 301 302 END SUBROUTINE hpg_zco 303 304 305 SUBROUTINE hpg_zps( kt ) 306 !!--------------------------------------------------------------------- 307 !! *** ROUTINE hpg_zps *** 308 !! 309 !! ** Method : z-coordinate plus partial steps case. blahblah... 310 !! 311 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 312 !! 313 !! History : 314 !! 8.5 ! 02-08 (A. Bozec) Original code 315 !! 9.0 ! 04-08 (G. Madec) F90 316 !!---------------------------------------------------------------------- 317 !! * modules used 318 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 319 & zhpj => sa ! use sa as 3D workspace 320 321 !! * Arguments 322 INTEGER, INTENT( in ) :: kt ! ocean time-step index 323 324 !! * local declarations 325 INTEGER :: ji, jj, jk ! dummy loop indices 326 INTEGER :: iku, ikv ! temporary integers 327 REAL(wp) :: & 328 zcoef0, zcoef1, & ! temporary scalars 329 zcoef2, zcoef3 ! " " 330 !!---------------------------------------------------------------------- 331 332 IF( kt == nit000 ) THEN 333 IF(lwp) WRITE(numout,*) 334 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 335 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps' 336 IF(lwp) WRITE(numout,*) ' vector optimization' 337 ENDIF 338 339 340 ! 0. Local constant initialization 341 ! -------------------------------- 342 zcoef0 = - grav * 0.5 343 344 ! 1. Surface value 345 ! ---------------- 346 DO jj = 2, jpjm1 347 DO ji = fs_2, fs_jpim1 ! vector opt. 348 zcoef1 = zcoef0 * fse3w(ji,jj,1) 349 ! hydrostatic pressure gradient 350 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 351 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 352 ! add to the general momentum trend 353 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 354 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 355 END DO 356 END DO 357 246 358 ! 2. interior value (2=<jk=<jpkm1) 247 359 ! ----------------- … … 252 364 ! hydrostatic pressure gradient 253 365 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 254 & + zcoef1 * ( ( rhd(ji+1,jj,jk) +rhd(ji+1,jj,jk-1) ) &255 & - ( rhd(ji ,jj,jk) +rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)366 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 367 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj) 256 368 257 369 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 258 & + zcoef1 * ( ( rhd(ji,jj+1,jk) +rhd(ji,jj+1,jk-1) ) &259 & - ( rhd(ji,jj, jk) +rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)370 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 371 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj) 260 372 ! add to the general momentum trend 261 373 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 262 374 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 263 END DO 375 END DO 264 376 END DO 265 377 END DO … … 279 391 ! on i-direction 280 392 IF ( iku > 2 ) THEN 281 ! subtract old value 393 ! subtract old value 282 394 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) 283 ! compute the new one 395 ! compute the new one 284 396 zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1) & 285 397 + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) … … 289 401 ! on j-direction 290 402 IF ( ikv > 2 ) THEN 291 ! subtract old value 403 ! subtract old value 292 404 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) 293 ! compute the new one 405 ! compute the new one 294 406 zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1) & 295 407 + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) … … 302 414 END DO 303 415 304 ! save the hydrostatic pressure gradient trends for diagnostic 305 ! momentum trends 306 IF( l_trddyn ) THEN 307 zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 308 zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 309 CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 310 ENDIF 311 312 IF(ln_ctl) THEN ! print sum trends (used for debugging) 313 CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, & 314 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 315 ENDIF 316 317 END SUBROUTINE dyn_hpg 318 319 #else 320 !!--------------------------------------------------------------------- 321 !! Default case : z-coordinate 322 !!--------------------------------------------------------------------- 323 324 SUBROUTINE dyn_hpg( kt ) 416 END SUBROUTINE hpg_zps 417 418 419 SUBROUTINE hpg_sco( kt ) 325 420 !!--------------------------------------------------------------------- 326 !! *** ROUTINE dyn_hpg *** 327 !! 328 !! ** Purpose : Compute the now momentum trend due to the horizontal 329 !! gradient of the hydrostatic pressure. Add it to the general 330 !! momentum trend. 331 !! 332 !! ** Method : The now hydrostatic pressure gradient at a given level 333 !! jk is computed by taking the vertical integral of the in-situ 421 !! *** ROUTINE hpg_sco *** 422 !! 423 !! ** Method : s-coordinate case. Jacobian scheme. 424 !! The now hydrostatic pressure gradient at a given level, jk, 425 !! is computed by taking the vertical integral of the in-situ 334 426 !! density gradient along the model level from the suface to that 335 !! level: zhpi = grav ..... 336 !! zhpj = grav ..... 427 !! level. s-coordinates (ln_sco): a corrective term is added 428 !! to the horizontal pressure gradient : 429 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 430 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 337 431 !! add it to the general momentum trend (ua,va). 338 !! 339 !! 432 !! ua = ua - 1/e1u * zhpi 433 !! va = va - 1/e2v * zhpj 340 434 !! 341 435 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 342 !! - Save the trend in (utrd,vtrd) ('key_trddyn')343 436 !! 344 437 !! History : 345 !! 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code 346 !! ! 91-11 (G. Madec) 347 !! ! 96-01 (G. Madec) s-coordinates 438 !! 7.0 ! 96-01 (G. Madec) s-coordinates 348 439 !! ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg 349 !! 8.5 ! 02-07 (G. Madec) F90: Free form and module 440 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module, vector opt. 441 !! 9.0 ! 04-08 (C. Talandier) New trends organization 442 !! 9.0 ! 05-10 (A. Beckmann) various s-coordinate options 350 443 !!---------------------------------------------------------------------- 351 444 !! * modules used … … 356 449 INTEGER, INTENT( in ) :: kt ! ocean time-step index 357 450 358 !! * local declarations451 !! * Local declarations 359 452 INTEGER :: ji, jj, jk ! dummy loop indices 360 453 REAL(wp) :: & 361 zcoef0, zcoef1, zuap, zvap ! temporary scalars 362 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 363 ztdua, ztdva ! temporary scalars 454 zcoef0, zuap, zvap ! temporary scalars 364 455 !!---------------------------------------------------------------------- 365 456 366 457 IF( kt == nit000 ) THEN 367 458 IF(lwp) WRITE(numout,*) 368 IF(lwp) WRITE(numout,*) 'dyn _hpg: hydrostatic pressure gradient trend'369 IF(lwp) WRITE(numout,*) '~~~~~~~ z-coordinate case'459 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 460 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 370 461 ENDIF 371 462 372 ! Save ua and va trends373 IF( l_trddyn ) THEN374 ztdua(:,:,:) = ua(:,:,:)375 ztdva(:,:,:) = va(:,:,:)376 ENDIF377 463 378 464 ! 0. Local constant initialization 379 465 ! -------------------------------- 380 466 zcoef0 = - grav * 0.5 381 zuap = 0.e0 382 zvap = 0.e0 383 467 468 469 ! 1. Surface value 470 ! ---------------- 471 DO jj = 2, jpjm1 472 DO ji = fs_2, fs_jpim1 ! vector opt. 473 ! hydrostatic pressure gradient along s-surfaces 474 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) & 475 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) ) 476 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) & 477 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) ) 478 ! s-coordinate pressure gradient correction 479 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 480 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 481 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 482 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 483 ! add to the general momentum trend 484 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 485 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 486 END DO 487 END DO 488 489 490 ! 2. interior value (2=<jk=<jpkm1) 491 ! ----------------- 492 DO jk = 2, jpkm1 493 DO jj = 2, jpjm1 494 DO ji = fs_2, fs_jpim1 ! vector opt. 495 ! hydrostatic pressure gradient along s-surfaces 496 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 497 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 498 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) 499 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 500 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 501 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) 502 ! s-coordinate pressure gradient correction 503 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 504 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 505 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 506 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 507 ! add to the general momentum trend 508 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 509 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 510 END DO 511 END DO 512 END DO 513 514 END SUBROUTINE hpg_sco 515 516 517 SUBROUTINE hpg_hel( kt ) 518 !!--------------------------------------------------------------------- 519 !! *** ROUTINE hpg_hel *** 520 !! 521 !! ** Method : s-coordinate case. 522 !! The now hydrostatic pressure gradient at a given level 523 !! jk is computed by taking the vertical integral of the in-situ 524 !! density gradient along the model level from the suface to that 525 !! level. s-coordinates (ln_sco): a corrective term is added 526 !! to the horizontal pressure gradient : 527 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 528 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 529 !! add it to the general momentum trend (ua,va). 530 !! ua = ua - 1/e1u * zhpi 531 !! va = va - 1/e2v * zhpj 532 !! 533 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 534 !! - Save the trend (l_trddyn=T) 535 !! 536 !! History : 537 !! 9.0 ! 05-10 (A. Beckmann) Original code 538 !!---------------------------------------------------------------------- 539 !! * modules used 540 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 541 & zhpj => sa ! use sa as 3D workspace 542 543 !! * Arguments 544 INTEGER, INTENT( in ) :: kt ! ocean time-step index 545 546 !! * Local declarations 547 INTEGER :: ji, jj, jk ! dummy loop indices 548 REAL(wp) :: & 549 zcoef0, zuap, zvap ! temporary scalars 550 !!---------------------------------------------------------------------- 551 552 IF( kt == nit000 ) THEN 553 IF(lwp) WRITE(numout,*) 554 IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend' 555 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, helsinki modified scheme' 556 ENDIF 557 558 559 ! 0. Local constant initialization 560 ! -------------------------------- 561 zcoef0 = - grav * 0.5 562 563 384 564 ! 1. Surface value 385 565 ! ---------------- 386 566 DO jj = 2, jpjm1 387 567 DO ji = fs_2, fs_jpim1 ! vector opt. 388 zcoef1 = zcoef0 * fse3w(ji,jj,1) 389 ! hydrostatic pressure gradient 390 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 391 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 568 ! hydrostatic pressure gradient along s-surfaces 569 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) & 570 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) ) 571 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) & 572 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) ) 573 ! s-coordinate pressure gradient correction 574 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 575 & * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 576 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 577 & * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 392 578 ! add to the general momentum trend 393 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 394 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 579 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 580 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 395 581 END DO 396 582 END DO … … 401 587 DO jj = 2, jpjm1 402 588 DO ji = fs_2, fs_jpim1 ! vector opt. 403 zcoef1 = zcoef0 * fse3w(ji,jj,jk) 404 ! hydrostatic pressure gradient 405 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 406 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 407 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj) 408 409 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 410 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 411 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj) 589 ! hydrostatic pressure gradient along s-surfaces 590 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 591 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk) & 592 & -fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk) ) & 593 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) & 594 & -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) ) 595 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 596 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) & 597 & -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) & 598 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 599 & -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) ) 600 ! s-coordinate pressure gradient correction 601 zuap = - zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 602 & * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 603 zvap = - zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 604 & * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 605 ! add to the general momentum trend 606 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 607 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 608 END DO 609 END DO 610 END DO 611 612 END SUBROUTINE hpg_hel 613 614 615 SUBROUTINE hpg_wdj( kt ) 616 !!--------------------------------------------------------------------- 617 !! *** ROUTINE hpg_wdj *** 618 !! 619 !! ** Method : Weighted Density Jacobian (wdj) scheme (song 1998) 620 !! The weighting coefficients from the namelist parameter gamm 621 !! (alpha=0.5-gamm ; beta=1-alpha=0.5+gamm) 622 !! 623 !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998. 624 !! 625 !! History : 626 !! 9.0 ! 05-05 (B.W. An) Original code 627 !! ! 05-10 (G. Madec) style & small optimisation 628 !!---------------------------------------------------------------------- 629 !! * modules used 630 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 631 & zhpj => sa ! use sa as 3D workspace 632 633 !! * Arguments 634 INTEGER, INTENT( in ) :: kt ! ocean time-step index 635 636 !! * Local declarations 637 INTEGER :: ji, jj, jk ! dummy loop indices 638 REAL(wp) :: & 639 zcoef0, zuap, zvap, & ! temporary scalars 640 zalph , zbeta ! " " 641 !!---------------------------------------------------------------------- 642 643 IF( kt == nit000 ) THEN 644 IF(lwp) WRITE(numout,*) 645 IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend' 646 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ Weighted Density Jacobian' 647 ENDIF 648 649 650 ! Local constant initialization 651 ! ----------------------------- 652 zcoef0 = - grav * 0.5 653 zalph = 0.5 - gamm ! weighting coefficients (alpha=0.5-gamm) 654 zbeta = 0.5 + gamm ! (beta =1-alpha=0.5+gamm) 655 656 ! Surface value (no ponderation) 657 ! ------------- 658 DO jj = 2, jpjm1 659 DO ji = fs_2, fs_jpim1 ! vector opt. 660 ! hydrostatic pressure gradient along s-surfaces 661 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) & 662 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) ) 663 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) & 664 & - fse3w(ji ,jj ,1) * rhd(ji, jj ,1) ) 665 ! s-coordinate pressure gradient correction 666 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 667 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 668 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 669 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 670 ! add to the general momentum trend 671 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 672 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 673 END DO 674 END DO 675 676 ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta) 677 ! -------------- 678 DO jk = 2, jpkm1 679 DO jj = 2, jpjm1 680 DO ji = fs_2, fs_jpim1 ! vector opt. 681 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 682 & * ( ( fsde3w(ji+1,jj,jk ) + fsde3w(ji,jj,jk ) & 683 & - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) & 684 & * ( zalph * ( rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) & 685 & + zbeta * ( rhd (ji+1,jj,jk ) - rhd (ji,jj,jk ) ) ) & 686 & - ( rhd (ji+1,jj,jk ) + rhd (ji,jj,jk ) & 687 & - rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) & 688 & * ( zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) & 689 & + zbeta * ( fsde3w(ji+1,jj,jk ) - fsde3w(ji,jj,jk ) ) ) ) 690 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 691 & * ( ( fsde3w(ji,jj+1,jk ) + fsde3w(ji,jj,jk ) & 692 & - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) & 693 & * ( zalph * ( rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) & 694 & + zbeta * ( rhd (ji,jj+1,jk ) - rhd (ji,jj,jk ) ) ) & 695 & - ( rhd (ji,jj+1,jk ) + rhd (ji,jj,jk ) & 696 & - rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) & 697 & * ( zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) & 698 & + zbeta * ( fsde3w(ji,jj+1,jk ) - fsde3w(ji,jj,jk ) ) ) ) 412 699 ! add to the general momentum trend 413 700 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 414 701 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 415 END DO 416 END DO 417 END DO 418 419 ! save the hydrostatic pressure ggradient trends for diagnostic 420 ! momentum trends 421 IF( l_trddyn ) THEN 422 zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 423 zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 424 425 CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 702 END DO 703 END DO 704 END DO 705 706 END SUBROUTINE hpg_wdj 707 708 709 SUBROUTINE hpg_djc( kt ) 710 !!--------------------------------------------------------------------- 711 !! *** ROUTINE hpg_djc *** 712 !! 713 !! ** Method : Density Jacobian with Cubic polynomial scheme 714 !! 715 !! Reference: Shchepetkin, A.F. & J.C. McWilliams, J. Geophys. Res., 716 !! 108(C3), 3090, 2003 717 !! History : 718 !! 9.0 ! 05-05 (B.W. An) Original code 719 !!---------------------------------------------------------------------- 720 !! * modules used 721 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 722 & zhpj => sa ! use sa as 3D workspace 723 724 !! * Arguments 725 INTEGER, INTENT( in ) :: kt ! ocean time-step index 726 727 !! * Local declarations 728 INTEGER :: ji, jj, jk ! dummy loop indices 729 REAL(wp) :: & 730 zcoef0, z1_10, cffu, cffx, & ! temporary scalars 731 z1_12, cffv, cffy, & ! " " 732 zep , cffw ! " " 733 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace 734 drhox, dzx, drhou, dzu, rho_i, & 735 drhoy, dzy, drhov, dzv, rho_j, & 736 drhoz, dzz, drhow, dzw, rho_k 737 !!---------------------------------------------------------------------- 738 739 IF( kt == nit000 ) THEN 740 IF(lwp) WRITE(numout,*) 741 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 742 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 426 743 ENDIF 427 744 428 IF(ln_ctl) THEN ! print sum trends (used for debugging) 429 CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, & 430 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 745 746 ! 0. Local constant initialization 747 ! -------------------------------- 748 zcoef0 = - grav * 0.5 749 z1_10 = 1.0 / 10.0 750 z1_12 = 1.0 / 12.0 751 752 !---------------------------------------------------------------------------------------- 753 ! compute and store in provisional arrays elementary vertical and horizontal differences 754 !---------------------------------------------------------------------------------------- 755 756 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 757 758 DO jk = 2, jpkm1 759 DO jj = 2, jpjm1 760 DO ji = fs_2, fs_jpim1 ! vector opt. 761 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 762 dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1) 763 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 764 dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk ) 765 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 766 dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk ) 767 END DO 768 END DO 769 END DO 770 771 !------------------------------------------------------------------------- 772 ! compute harmonic averages using eq. 5.18 773 !------------------------------------------------------------------------- 774 zep = 1.e-15 775 776 !!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2) 777 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 778 779 DO jk = 2, jpkm1 780 DO jj = 2, jpjm1 781 DO ji = fs_2, fs_jpim1 ! vector opt. 782 cffw = 2.0 * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 783 784 cffu = 2.0 * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 785 cffx = 2.0 * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 786 787 cffv = 2.0 * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 788 cffy = 2.0 * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 789 790 IF( cffw > zep) THEN 791 drhow(ji,jj,jk) = 2.0 * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 792 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 793 ELSE 794 drhow(ji,jj,jk) = 0.e0 795 ENDIF 796 797 dzw(ji,jj,jk) = 2.0 * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 798 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 799 800 IF( cffu > zep ) THEN 801 drhou(ji,jj,jk) = 2.0 * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 802 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 803 ELSE 804 drhou(ji,jj,jk ) = 0.e0 805 ENDIF 806 807 IF( cffx > zep ) THEN 808 dzu(ji,jj,jk) = 2.0*dzx(ji+1,jj,jk)*dzx(ji,jj,jk) & 809 & /(dzx(ji+1,jj,jk)+dzx(ji,jj,jk)) 810 ELSE 811 dzu(ji,jj,jk) = 0.e0 812 ENDIF 813 814 IF( cffv > zep ) THEN 815 drhov(ji,jj,jk) = 2.0 * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 816 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 817 ELSE 818 drhov(ji,jj,jk) = 0.e0 819 ENDIF 820 821 IF( cffy > zep ) THEN 822 dzv(ji,jj,jk) = 2.0 * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 823 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 824 ELSE 825 dzv(ji,jj,jk) = 0.e0 826 ENDIF 827 828 END DO 829 END DO 830 END DO 831 832 !---------------------------------------------------------------------------------- 833 ! apply boundary conditions at top and bottom using 5.36-5.37 834 !---------------------------------------------------------------------------------- 835 drhow(:,:, 1 ) = 1.5 * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5 * drhow(:,:, 2 ) 836 drhou(:,:, 1 ) = 1.5 * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5 * drhou(:,:, 2 ) 837 drhov(:,:, 1 ) = 1.5 * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5 * drhov(:,:, 2 ) 838 839 drhow(:,:,jpk) = 1.5 * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5 * drhow(:,:,jpkm1) 840 drhou(:,:,jpk) = 1.5 * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5 * drhou(:,:,jpkm1) 841 drhov(:,:,jpk) = 1.5 * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5 * drhov(:,:,jpkm1) 842 843 844 !-------------------------------------------------------------- 845 ! Upper half of top-most grid box, compute and store 846 !------------------------------------------------------------- 847 848 !!bug gm : e3w-de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified 849 ! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 850 851 DO jj = 2, jpjm1 852 DO ji = fs_2, fs_jpim1 ! vector opt. 853 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) & 854 & * ( rhd(ji,jj,1) & 855 & + 0.5 * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) & 856 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) & 857 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) ) 858 END DO 859 END DO 860 861 !!bug gm : here also, simplification is possible 862 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop 863 864 DO jk = 2, jpkm1 865 DO jj = 2, jpjm1 866 DO ji = fs_2, fs_jpim1 ! vector opt. 867 868 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 869 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) & 870 & - grav * z1_10 * ( & 871 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 872 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 873 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 874 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 875 & ) 876 877 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 878 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) & 879 & - grav* z1_10 * ( & 880 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 881 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 882 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 883 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 884 & ) 885 886 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 887 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) & 888 & - grav* z1_10 * ( & 889 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 890 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 891 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 892 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 893 & ) 894 895 END DO 896 END DO 897 END DO 898 CALL lbc_lnk(rho_k,'W',1.) 899 CALL lbc_lnk(rho_i,'U',1.) 900 CALL lbc_lnk(rho_j,'V',1.) 901 902 903 ! --------------- 904 ! Surface value 905 ! --------------- 906 DO jj = 2, jpjm1 907 DO ji = fs_2, fs_jpim1 ! vector opt. 908 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 909 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 910 ! add to the general momentum trend 911 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 912 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 913 END DO 914 END DO 915 916 ! ---------------- 917 ! interior value (2=<jk=<jpkm1) 918 ! ---------------- 919 DO jk = 2, jpkm1 920 DO jj = 2, jpjm1 921 DO ji = fs_2, fs_jpim1 ! vector opt. 922 ! hydrostatic pressure gradient along s-surfaces 923 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 924 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 925 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) / e1u(ji,jj) 926 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 927 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 928 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 929 ! add to the general momentum trend 930 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 931 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 932 END DO 933 END DO 934 END DO 935 936 END SUBROUTINE hpg_djc 937 938 939 SUBROUTINE hpg_rot( kt ) 940 !!--------------------------------------------------------------------- 941 !! *** ROUTINE hpg_rot *** 942 !! 943 !! ** Method : rotated axes scheme (Thiem and Berntsen 2005) 944 !! 945 !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005. 946 !! History : 947 !! 9.0 ! 05-07 (B.W. An) 948 !! 9.0 ! 05-10 (A. Beckmann) adapted to non-equidistant and masked grids 949 !!---------------------------------------------------------------------- 950 !! * modules used 951 USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace 952 & zhpj => sa ! use sa as 3D workspace 953 954 !! * Arguments 955 INTEGER, INTENT( in ) :: kt ! ocean time-step index 956 957 !! * Local declarations 958 INTEGER :: ji, jj, jk ! dummy loop indices 959 REAL(wp) :: & 960 zforg, zcoef0, zuap, zmskd1, zmskd1m, & 961 zfrot , zvap, zmskd2, zmskd2m 962 REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D temporary workspace 963 zdistr, zsina, zcosa 964 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D temporary workspace 965 zhpiorg, zhpirot, zhpitra, zhpine, & 966 zhpjorg, zhpjrot, zhpjtra, zhpjne 967 !!---------------------------------------------------------------------- 968 969 IF( kt == nit000 ) THEN 970 IF(lwp) WRITE(numout,*) 971 IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend' 972 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, rotated axes scheme used' 431 973 ENDIF 432 974 433 END SUBROUTINE dyn_hpg 434 435 #endif 975 ! ------------------------------- 976 ! Local constant initialization 977 ! ------------------------------- 978 zcoef0 = - grav * 0.5 979 zforg = 0.95e0 980 zfrot = 1.e0 - zforg 981 982 ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 983 zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) ) 984 985 ! sinus and cosinus of diagonal angle at F-point 986 zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) ) 987 zcosa(:,:) = COS( zsina(:,:) ) 988 zsina(:,:) = SIN( zsina(:,:) ) 989 990 ! --------------- 991 ! Surface value 992 ! --------------- 993 ! compute and add to the general trend the pressure gradients along the axes 994 DO jj = 2, jpjm1 995 DO ji = fs_2, fs_jpim1 ! vector opt. 996 ! hydrostatic pressure gradient along s-surfaces 997 zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,1) * rhd(ji+1,jj,1) & 998 & - fse3t(ji ,jj,1) * rhd(ji ,jj,1) ) 999 zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,1) * rhd(ji,jj+1,1) & 1000 & - fse3t(ji,jj ,1) * rhd(ji,jj ,1) ) 1001 ! s-coordinate pressure gradient correction 1002 zuap = -zcoef0 * ( rhd (ji+1,jj ,1) + rhd (ji,jj,1) ) & 1003 & * ( fsdept(ji+1,jj ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 1004 zvap = -zcoef0 * ( rhd (ji ,jj+1,1) + rhd (ji,jj,1) ) & 1005 & * ( fsdept(ji ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 1006 ! add to the general momentum trend 1007 ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap ) 1008 va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap ) 1009 END DO 1010 END DO 1011 1012 ! compute the pressure gradients in the diagonal directions 1013 DO jj = 1, jpjm1 1014 DO ji = 1, fs_jpim1 ! vector opt. 1015 zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji ,jj,1) ! mask in the 1st diagnonal 1016 zmskd2 = tmask(ji ,jj+1,1) * tmask(ji+1,jj,1) ! mask in the 2nd diagnonal 1017 ! hydrostatic pressure gradient along s-surfaces 1018 zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1) & 1019 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) ) 1020 zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) & 1021 & - fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) ) 1022 ! s-coordinate pressure gradient correction 1023 zuap = -zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,1) + rhd (ji ,jj,1) ) & 1024 & * ( fsdept(ji+1,jj+1,1) - fsdept(ji ,jj,1) ) 1025 zvap = -zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,1) + rhd (ji+1,jj,1) ) & 1026 & * ( fsdept(ji ,jj+1,1) - fsdept(ji+1,jj,1) ) 1027 ! back rotation 1028 zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) & 1029 & - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 1030 zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) & 1031 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 1032 END DO 1033 END DO 1034 1035 ! interpolate and add to the general trend the diagonal gradient 1036 DO jj = 2, jpjm1 1037 DO ji = fs_2, fs_jpim1 ! vector opt. 1038 ! averaging 1039 zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji ,jj-1,1) ) 1040 zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj ,1) ) 1041 ! add to the general momentum trend 1042 ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1) 1043 va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1) 1044 END DO 1045 END DO 1046 1047 ! ----------------- 1048 ! 2. interior value (2=<jk=<jpkm1) 1049 ! ----------------- 1050 ! compute and add to the general trend the pressure gradients along the axes 1051 DO jk = 2, jpkm1 1052 DO jj = 2, jpjm1 1053 DO ji = fs_2, fs_jpim1 ! vector opt. 1054 ! hydrostatic pressure gradient along s-surfaces 1055 zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1) & 1056 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk ) & 1057 & - fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk ) & 1058 & + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) & 1059 & - fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) ) 1060 zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1) & 1061 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk ) & 1062 & - fse3t(ji,jj ,jk ) * rhd(ji,jj, jk ) & 1063 & + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 1064 & - fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) ) 1065 ! s-coordinate pressure gradient correction 1066 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 1067 & * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 1068 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 1069 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 1070 ! add to the general momentum trend 1071 ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap ) 1072 va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap ) 1073 END DO 1074 END DO 1075 END DO 1076 1077 ! compute the pressure gradients in the diagonal directions 1078 DO jk = 2, jpkm1 1079 DO jj = 1, jpjm1 1080 DO ji = 1, fs_jpim1 ! vector opt. 1081 zmskd1 = tmask(ji+1,jj+1,jk ) * tmask(ji ,jj,jk ) ! level jk mask in the 1st diagnonal 1082 zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji ,jj,jk-1) ! level jk-1 " " 1083 zmskd2 = tmask(ji ,jj+1,jk ) * tmask(ji+1,jj,jk ) ! level jk mask in the 2nd diagnonal 1084 zmskd2m = tmask(ji ,jj+1,jk-1) * tmask(ji+1,jj,jk-1) ! level jk-1 " " 1085 ! hydrostatic pressure gradient along s-surfaces 1086 zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1) & 1087 & + zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,jk ) * rhd(ji+1,jj+1,jk) & 1088 & -fse3t(ji ,jj ,jk ) * rhd(ji ,jj ,jk) ) & 1089 & + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1) & 1090 & -fse3t(ji ,jj ,jk-1) * rhd(ji ,jj ,jk-1) ) 1091 zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1) & 1092 & + zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,jk ) * rhd(ji ,jj+1,jk) & 1093 & -fse3t(ji+1,jj ,jk ) * rhd(ji+1,jj, jk) ) & 1094 & + zdistr(ji,jj) * zmskd2m * ( fse3t(ji ,jj+1,jk-1) * rhd(ji ,jj+1,jk-1) & 1095 & -fse3t(ji+1,jj ,jk-1) * rhd(ji+1,jj, jk-1) ) 1096 ! s-coordinate pressure gradient correction 1097 zuap = - zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,jk) + rhd (ji ,jj,jk) ) & 1098 & * ( fsdept(ji+1,jj+1,jk) - fsdept(ji ,jj,jk) ) 1099 zvap = - zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,jk) + rhd (ji+1,jj,jk) ) & 1100 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji+1,jj,jk) ) 1101 ! back rotation 1102 zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) & 1103 & - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 1104 zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) & 1105 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 1106 END DO 1107 END DO 1108 END DO 1109 1110 ! interpolate and add to the general trend 1111 DO jk = 2, jpkm1 1112 DO jj = 2, jpjm1 1113 DO ji = fs_2, fs_jpim1 ! vector opt. 1114 ! averaging 1115 zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji ,jj-1,jk) ) 1116 zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj ,jk) ) 1117 ! add to the general momentum trend 1118 ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk) 1119 va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk) 1120 END DO 1121 END DO 1122 END DO 1123 1124 END SUBROUTINE hpg_rot 436 1125 437 1126 !!====================================================================== -
trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r258 r455 17 17 USE trdmod_oce ! ocean variables trends 18 18 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 USE prtctl ! Print control20 19 21 20 IMPLICIT NONE … … 61 60 !! diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ] 62 61 !! diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ] 63 !! If l k_sco=F and lk_zps=F, the vertical scale factors in the62 !! If ln_sco=F and ln_zps=F, the vertical scale factors in the 64 63 !! rotational part of the diffusion are simplified 65 64 !! Add this before trend to the general trend (ua,va): … … 70 69 !! ** Action : - Update (ua,va) with the before iso-level biharmonic 71 70 !! mixing trend. 72 !! - Save in (ztdua,ztdva) the trends ('key_trddyn')73 71 !! 74 72 !! History : … … 81 79 !! 9.0 ! 04-08 (C. Talandier) New trends organization 82 80 !!---------------------------------------------------------------------- 83 !! * Modules used84 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace85 ztdva => sa ! use sa as 3D workspace86 87 81 !! * Arguments 88 82 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 107 101 zlv(:,:) = 0.e0 108 102 109 ! Save ua and va trends110 IF( l_trddyn ) THEN111 ztdua(:,:,:) = ua(:,:,:)112 ztdva(:,:,:) = va(:,:,:)113 ENDIF114 103 ! ! =============== 115 104 DO jk = 1, jpkm1 ! Horizontal slab … … 118 107 ! --------- 119 108 120 IF( l k_sco .OR. lk_zps ) THEN ! s-coordinate or z-coordinate with partial steps109 IF( ln_sco .OR. ln_zps ) THEN ! s-coordinate or z-coordinate with partial steps 121 110 zuf(:,:) = rotb(:,:,jk) * fse3f(:,:,jk) 122 111 DO jj = 2, jpjm1 … … 129 118 END DO 130 119 END DO 131 ELSE ! z-coordinate 120 ELSE ! z-coordinate - full step 132 121 DO jj = 2, jpjm1 133 122 DO ji = fs_2, fs_jpim1 ! vector opt. … … 162 151 zuf(ji,jj) = fmask(ji,jj,jk) * ( zcv(ji+1,jj ) - zcv(ji,jj) & 163 152 & - zcu(ji ,jj+1) + zcu(ji,jj) ) & 164 #if defined key_s_coord || defined key_partial_steps 153 #if defined key_zco 154 & / ( e1f(ji,jj)*e2f(ji,jj) ) 155 #else 165 156 & * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) 166 #else167 & / ( e1f(ji,jj)*e2f(ji,jj) )168 157 #endif 169 158 END DO … … 173 162 DO jj = 1, jpjm1 174 163 DO ji = 1, fs_jpim1 ! vector opt. 175 #if defined key_s_coord || defined key_partial_steps 164 #if defined key_zco 165 zlu(ji,jj) = e2u(ji,jj) * zlu(ji,jj) 166 zlv(ji,jj) = e1v(ji,jj) * zlv(ji,jj) 167 #else 176 168 zlu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj) 177 169 zlv(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj) 178 #else179 zlu(ji,jj) = e2u(ji,jj) * zlu(ji,jj)180 zlv(ji,jj) = e1v(ji,jj) * zlv(ji,jj)181 170 #endif 182 171 END DO … … 186 175 DO jj = 2, jpj 187 176 DO ji = fs_2, jpi ! vector opt. 188 #if defined key_s_coord || defined key_partial_steps 177 #if defined key_zco 178 zbt = e1t(ji,jj) * e2t(ji,jj) 179 #else 189 180 zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 190 #else191 zbt = e1t(ji,jj) * e2t(ji,jj)192 181 #endif 193 182 zut(ji,jj) = ( zlu(ji,jj) - zlu(ji-1,jj ) & … … 207 196 DO jj = 2, jpjm1 208 197 DO ji = fs_2, fs_jpim1 ! vector opt. 209 #if defined key_s_coord || defined key_partial_steps 198 #if defined key_zco 199 ze2u = e2u(ji,jj) 200 ze2v = e1v(ji,jj) 201 #else 210 202 ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) 211 203 ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) 212 #else213 ze2u = e2u(ji,jj)214 ze2v = e1v(ji,jj)215 204 #endif 216 205 ! horizontal biharmonic diffusive trends … … 229 218 END DO ! End of slab 230 219 ! ! =============== 231 ! save the lateral diffusion trends for diagnostic232 ! momentum trends233 IF( l_trddyn ) THEN234 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)235 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)236 237 CALL trd_mod(ztdua, ztdva, jpdtdldf, 'DYN', kt)238 ENDIF239 240 IF(ln_ctl) THEN ! print sum trends (used for debugging)241 CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf - Ua: ', mask1=umask, &242 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')243 ENDIF244 220 245 221 END SUBROUTINE dyn_ldf_bilap -
trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r258 r455 130 130 END DO ! End of slab 131 131 ! ! =============== 132 133 ! save the lateral diffusion trends for diagnostic134 ! momentum trends135 IF( l_trddyn ) THEN136 137 CALL trd_mod(zwk3, zwk4, jpdtdldf, 'DYN', kt)138 ENDIF139 140 IF(ln_ctl) THEN ! print sum trends (used for debugging)141 CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf - Ua: ', mask1=umask, &142 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')143 ENDIF144 132 145 133 END SUBROUTINE dyn_ldf_bilapg -
trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r258 r455 21 21 USE trdmod_oce ! ocean variables trends 22 22 USE ldfslp ! iso-neutral slopes 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 24 USE in_out_manager ! I/O manager 24 25 USE prtctl ! Print control … … 36 37 !!---------------------------------------------------------------------- 37 38 !! OPA 9.0 , LOCEAN-IPSL (2005) 38 !! $Header$39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt40 39 !!---------------------------------------------------------------------- 41 40 … … 46 45 !! *** ROUTINE dyn_ldf_iso *** 47 46 !! 48 !! ** Purpose : Compute the before trend of the horizontal part of the 49 !! lateral momentum diffusion and add it to the general trend of 50 !! momentum equation. 47 !! ** Purpose : Compute the before trend of the rotated laplacian 48 !! operator of lateral momentum diffusion except the diagonal 49 !! vertical term that will be computed in dynzdf module. Add it 50 !! to the general trend of momentum equation. 51 51 !! 52 52 !! ** Method : 53 !! The horizontal component of the lateral diffusive trends on54 !! momentum is provided by a 2nd order operator rotated along neu-55 !! tral or geopotential surfaces(in s-coordinates).53 !! The momentum lateral diffusive trend is provided by a 2nd 54 !! order operator rotated along neutral or geopotential surfaces 55 !! (in s-coordinates). 56 56 !! It is computed using before fields (forward in time) and isopyc- 57 !! nal or geopotential slopes computed in routine ldfslp or inildf.57 !! nal or geopotential slopes computed in routine ldfslp. 58 58 !! Here, u and v components are considered as 2 independent scalar 59 59 !! fields. Therefore, the property of splitting divergent and rota- … … 76 76 !! Add this trend to the general trend (ua,va): 77 77 !! ua = ua + diffu 78 !! 'key_trddyn' defined: the trends are saved for diagnostics. 78 !! CAUTION: here the isopycnal part is with a coeff. of aht. This 79 !! should be modified for applications others than orca_r2 (!!bug) 79 80 !! 80 81 !! ** Action : 81 82 !! Update (ua,va) arrays with the before geopotential biharmonic 82 83 !! mixing trend. 83 !! Save in (uldftrd,vldftrd) arrays the trends if 'key_trddyn' defined 84 !! Update (avmu,avmv) to accompt for the diagonal vertical component 85 !! of the rotated operator in dynzdf module 84 86 !! 85 87 !! History : … … 87 89 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 88 90 !! 9.0 ! 04-08 (C. Talandier) New trends organization 91 !! ! 05-11 (G. Madec) s-coordinate: horizontal diffusion 89 92 !!---------------------------------------------------------------------- 90 !! * Modules used91 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace92 ztdva => sa ! use sa as 3D workspace93 93 !! * Arguments 94 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 103 103 ziut, zjuf, zjvt, zivf, & ! temporary workspace 104 104 zdku, zdk1u, zdkv, zdk1v 105 106 ! ajout dynzdf_iso 107 REAL(wp) :: & 108 zcoef0, zcoef3, zcoef4, zmkt, zmkf, & 109 zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj 110 REAL(wp), DIMENSION(jpi,jpj) :: & 111 zfuw, zdiu, zdju, zdj1u, & ! " " 112 zfvw, zdiv, zdjv, zdj1v 113 105 114 !!---------------------------------------------------------------------- 106 115 … … 111 120 ENDIF 112 121 113 ! Save ua and va trends 114 IF( l_trddyn ) THEN 115 ztdua(:,:,:) = ua(:,:,:) 116 ztdva(:,:,:) = va(:,:,:) 122 ! ! s-coordinate: Iso-level diffusion on momentum but not on tracer 123 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 124 125 ! set the slopes of iso-level 126 DO jk = 1, jpk 127 DO jj = 2, jpjm1 128 DO ji = fs_2, fs_jpim1 ! vector opt. 129 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 130 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 131 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 132 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 133 END DO 134 END DO 135 END DO 136 137 ! Lateral boundary conditions on the slopes 138 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 139 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 140 141 !!bug 142 if( kt == nit000 ) then 143 IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 144 & ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 145 endif 146 !!end 117 147 ENDIF 148 118 149 ! ! =============== 119 150 DO jk = 1, jpkm1 ! Horizontal slab … … 142 173 ! i-flux at t-point -----f----- 143 174 144 DO jj = 2, jpjm1 145 DO ji = fs_2, jpi ! vector opt. 146 zabe1 = ( fsahmt(ji,jj,jk) + ahmb0 ) & 147 #if defined key_partial_steps 148 * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1, jj,jk) ) / e1t(ji,jj) 149 #else 150 * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 151 #endif 152 153 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 154 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) 155 156 zcof1 = - aht0 * e2t(ji,jj) * zmskt & 157 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 158 159 ziut(ji,jj) = tmask(ji,jj,jk) * & 160 ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 161 + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 162 +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) 163 END DO 164 END DO 175 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 176 DO jj = 2, jpjm1 177 DO ji = fs_2, jpi ! vector opt. 178 zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 179 180 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 181 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) 182 183 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 184 185 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 186 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 187 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) 188 END DO 189 END DO 190 ELSE ! other coordinate system (zco or sco) : e3t 191 DO jj = 2, jpjm1 192 DO ji = fs_2, jpi ! vector opt. 193 zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 194 195 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 196 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) 197 198 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 199 200 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 201 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 202 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) 203 END DO 204 END DO 205 ENDIF 165 206 166 207 ! j-flux at f-point 167 208 DO jj = 1, jpjm1 168 209 DO ji = 1, fs_jpim1 ! vector opt. 169 zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) & 170 * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 210 zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 171 211 172 212 zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 173 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) 174 175 zcof2 = - aht0 * e1f(ji,jj) * zmskf & 176 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 177 178 zjuf(ji,jj) = fmask(ji,jj,jk) * & 179 ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & 180 + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 181 +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) 213 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) 214 215 zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 216 217 zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & 218 & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 219 & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk) 182 220 END DO 183 221 END DO … … 191 229 DO jj = 2, jpjm1 192 230 DO ji = 1, fs_jpim1 ! vector opt. 193 zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) & 194 * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 231 zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 195 232 196 233 zmskf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 197 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) 198 199 zcof1 = - aht0 * e2f(ji,jj) * zmskf & 200 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 201 202 zivf(ji,jj) = fmask(ji,jj,jk) * & 203 ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 204 + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 205 +zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) 234 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) 235 236 zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 237 238 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 239 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 240 & +zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) 206 241 END DO 207 242 END DO 208 243 209 244 ! j-flux at t-point 210 DO jj = 2, jpj 211 DO ji = 1, fs_jpim1 ! vector opt. 212 zabe2 = ( fsahmt(ji,jj,jk) + ahmb0 ) & 213 #if defined key_partial_steps 214 * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji, jj-1, jk) ) / e2t(ji,jj) 215 #else 216 * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 217 #endif 218 219 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 220 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 221 222 zcof2 = - aht0 * e1t(ji,jj) * zmskt & 223 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 224 225 zjvt(ji,jj) = tmask(ji,jj,jk) * & 226 ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 227 + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 228 +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) 229 END DO 230 END DO 245 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 246 DO jj = 2, jpj 247 DO ji = 1, fs_jpim1 ! vector opt. 248 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 249 250 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 251 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 252 253 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 254 255 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 256 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 257 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) 258 END DO 259 END DO 260 ELSE ! other coordinate system (zco or sco) : e3t 261 DO jj = 2, jpj 262 DO ji = 1, fs_jpim1 ! vector opt. 263 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 264 265 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 266 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 267 268 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 269 270 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 271 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 272 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) 273 END DO 274 END DO 275 ENDIF 231 276 232 277 … … 241 286 ! horizontal component of isopycnal momentum diffusive trends 242 287 zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + & 243 288 & zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu 244 289 zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + & 245 290 & zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv 246 291 ! add the trends to the general trends 247 292 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah … … 253 298 ! ! =============== 254 299 255 ! save the lateral diffusion trends for diagnostic 256 ! momentum trends will be saved in dynzdf_iso.F90 257 IF( l_trddyn ) THEN 258 uldftrd(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 259 vldftrd(:,:,:) = va(:,:,:) - ztdva(:,:,:) 260 ENDIF 261 262 IF(ln_ctl) THEN ! print sum trends (used for debugging) 263 CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf - Ua: ', mask1=umask, & 264 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 265 ENDIF 300 ! print sum trends (used for debugging) 301 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, & 302 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 303 304 305 ! ! =============== 306 DO jj = 2, jpjm1 ! Vertical slab 307 ! ! =============== 308 309 310 ! I. vertical trends associated with the lateral mixing 311 ! ===================================================== 312 ! (excluding the vertical flux proportional to dk[t] 313 314 315 ! I.1 horizontal momentum gradient 316 ! -------------------------------- 317 318 DO jk = 1, jpk 319 DO ji = 2, jpi 320 ! i-gradient of u at jj 321 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) ) 322 ! j-gradient of u and v at jj 323 zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) ) 324 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) ) 325 ! j-gradient of u and v at jj+1 326 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) ) 327 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) ) 328 END DO 329 END DO 330 DO jk = 1, jpk 331 DO ji = 1, jpim1 332 ! i-gradient of v at jj 333 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) ) 334 END DO 335 END DO 336 337 338 ! I.2 Vertical fluxes 339 ! ------------------- 340 341 ! Surface and bottom vertical fluxes set to zero 342 DO ji = 1, jpi 343 zfuw(ji, 1 ) = 0.e0 344 zfvw(ji, 1 ) = 0.e0 345 zfuw(ji,jpk) = 0.e0 346 zfvw(ji,jpk) = 0.e0 347 END DO 348 349 ! interior (2=<jk=<jpk-1) on U field 350 DO jk = 2, jpkm1 351 DO ji = 2, jpim1 352 zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 353 354 zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 355 zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 356 357 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 358 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. ) 359 zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) & 360 + fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. ) 361 362 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 363 zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj 364 ! vertical flux on u field 365 zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) & 366 +zdiu (ji,jk ) + zdiu (ji+1,jk ) ) & 367 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) & 368 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 369 ! update avmu (add isopycnal vertical coefficient to avmu) 370 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 371 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 372 END DO 373 END DO 374 375 ! interior (2=<jk=<jpk-1) on V field 376 DO jk = 2, jpkm1 377 DO ji = 2, jpim1 378 zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 379 380 zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 381 zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 382 383 zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) & 384 + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. ) 385 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) & 386 + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. ) 387 388 zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi 389 zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj 390 ! vertical flux on v field 391 zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) & 392 +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) & 393 + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 394 +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 395 ! update avmv (add isopycnal vertical coefficient to avmv) 396 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 397 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 398 END DO 399 END DO 400 401 402 ! I.3 Divergence of vertical fluxes added to the general tracer trend 403 ! ------------------------------------------------------------------- 404 DO jk = 1, jpkm1 405 DO ji = 2, jpim1 406 ! volume elements 407 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 408 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 409 ! part of the k-component of isopycnal momentum diffusive trends 410 zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 411 zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 412 ! add the trends to the general trends 413 ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 414 va(ji,jj,jk) = va(ji,jj,jk) + zvav 415 END DO 416 END DO 417 ! ! =============== 418 END DO ! End of slab 419 ! ! =============== 266 420 267 421 END SUBROUTINE dyn_ldf_iso -
trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90
r258 r455 18 18 USE trdmod_oce ! ocean variables trends 19 19 USE ldfslp ! iso-neutral slopes 20 USE prtctl ! Print control21 20 22 21 IMPLICIT NONE … … 51 50 !! difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb] 52 51 !! difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb] 53 !! If 'key_s_coord' key is not activated, the vertical scale factor54 !! i s simplified in the rotational part of the diffusion.52 !! If lk_zco=T, e3f=e3u=e3v, the vertical scale factor are simplified 53 !! in the rotational part of the diffusion. 55 54 !! Add this before trend to the general trend (ua,va): 56 55 !! (ua,va) = (ua,va) + (diffu,diffv) … … 60 59 !! ** Action : - Update (ua,va) with the before iso-level harmonic 61 60 !! mixing trend. 62 !! - Save in (ztdua,ztdva) arrays the trends ('key_trddyn')63 61 !! 64 62 !! History : … … 69 67 !! 9.0 ! 04-08 (C. Talandier) New trends organization 70 68 !!---------------------------------------------------------------------- 71 !! * Modules used72 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace73 ztdva => sa ! use sa as 3D workspace74 75 69 !! * Arguments 76 70 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 88 82 ENDIF 89 83 90 ! Save ua and va trends91 IF( l_trddyn ) THEN92 ztdua(:,:,:) = ua(:,:,:)93 ztdva(:,:,:) = va(:,:,:)94 ENDIF95 96 84 ! ! =============== 97 85 DO jk = 1, jpkm1 ! Horizontal slab … … 99 87 DO jj = 2, jpjm1 100 88 DO ji = fs_2, fs_jpim1 ! vector opt. 101 #if defined key_s_coord || defined key_partial_steps 89 #if defined key_zco 90 ! horizontal diffusive trends 91 ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk) 92 ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk) 93 zua = - ( ze2u - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk) ) / e2u(ji,jj) & 94 + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) - ze1v ) / e1u(ji,jj) 95 96 zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk) ) / e1v(ji,jj) & 97 + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v ) / e2v(ji,jj) 98 #else 102 99 ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk)*fse3f(ji,jj,jk) 103 100 ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk) … … 108 105 zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk)*fse3f(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) & 109 106 + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v ) / e2v(ji,jj) 110 #else111 ! horizontal diffusive trends112 ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk)113 ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk)114 zua = - ( ze2u - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk) ) / e2u(ji,jj) &115 + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) - ze1v ) / e1u(ji,jj)116 117 zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk) ) / e1v(ji,jj) &118 + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v ) / e2v(ji,jj)119 107 #endif 120 108 … … 128 116 ! ! =============== 129 117 130 ! save the lateral diffusion trends for diagnostic131 ! momentum trends132 IF( l_trddyn ) THEN133 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)134 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)135 136 CALL trd_mod(ztdua, ztdva, jpdtdldf, 'DYN', kt)137 ENDIF138 139 IF(ln_ctl) THEN ! print sum trends (used for debugging)140 CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf - Ua: ', mask1=umask, &141 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')142 ENDIF143 144 118 END SUBROUTINE dyn_ldf_lap 145 119 -
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r371 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_exp && ! defined key_ autotasking) || defined key_esopa7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_exp' free sfce cst vol. without filter nor ts9 !! NOT 'key_ autotasking'k-j-i loop (vector opt.)6 #if ( defined key_dynspg_exp && ! defined key_mpp_omp ) || defined key_esopa 7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_exp' free sfce cst vol. without filter nor ts 9 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_exp : update the momentum trend with the surface -
trunk/NEMO/OPA_SRC/DYN/dynspg_exp_jki.F90
r371 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_exp && defined key_ autotasking) || defined key_esopa6 #if ( defined key_dynspg_exp && defined key_mpp_omp ) || defined key_esopa 7 7 !!---------------------------------------------------------------------- 8 8 !! 'key_dynspg_exp' explicit free surface 9 !! 'key_ autotasking'j-k-i loop9 !! 'key_mpp_omp' j-k-i loop 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_exp_jki : update the momentum trend with the surface -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r413 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_flt && ! defined key_ autotasking) || defined key_esopa6 #if ( defined key_dynspg_flt && ! defined key_mpp_omp ) || defined key_esopa 7 7 !!---------------------------------------------------------------------- 8 8 !! 'key_dynspg_flt' filtered free surface 9 !! NOT 'key_ autotasking'k-j-i loop (vector opt.)9 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_flt : update the momentum trend with the surface pressure -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
r413 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_flt && defined key_ autotasking) || defined key_esopa7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_flt' & 'key_autotasking'filtered free surface9 !! j-k-i loop (j-slab)6 #if ( defined key_dynspg_flt && defined key_mpp_omp ) || defined key_esopa 7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_flt' filtered free surface 9 !! 'key_mpp_omp' j-k-i loop (vector opt.) 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_flt_jki : Update the momentum trend with the surface pressure -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r374 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_ts && ! defined key_ autotasking) || defined key_esopa7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_ts' free surface cst volume with time splitting9 !! NOT 'key_ autotasking'k-j-i loop (vector opt.)6 #if ( defined key_dynspg_ts && ! defined key_mpp_omp ) || defined key_esopa 7 !!---------------------------------------------------------------------- 8 !! 'key_dynspg_ts' free surface cst volume with time splitting 9 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_ts : compute surface pressure gradient trend using a time- -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts_jki.F90
r374 r455 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_ts && defined key_ autotasking) || defined key_esopa6 #if ( defined key_dynspg_ts && defined key_mpp_omp ) || defined key_esopa 7 7 !!---------------------------------------------------------------------- 8 8 !! 'key_dynspg_ts' free surface with time splitting 9 !! 'key_ autotasking'j-k-i loop (vector opt.)9 !! 'key_mpp_omp' j-k-i loop (vector opt.) 10 10 !!---------------------------------------------------------------------- 11 11 !! dyn_spg_ts : compute surface pressure gradient trend using a time- -
trunk/NEMO/OPA_SRC/DYN/dynvor.F90
r418 r455 7 7 8 8 !!---------------------------------------------------------------------- 9 !! dyn_vor_enstrophy: enstrophy conserving scheme (ln_dynvor_ens=T) 10 !! dyn_vor_energy : energy conserving scheme (ln_dynvor_ene=T) 11 !! dyn_vor_mixed : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 12 !! dyn_vor_ene_ens : energy and enstrophy conserving (ln_dynvor_een=T) 13 !! dyn_vor_ctl : control of the different vorticity option 9 !! dyn_vor : Update the momentum trend with the vorticity trend 10 !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) 11 !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) 12 !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 13 !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) 14 !! vor_ctl : control of the different vorticity option 14 15 !!---------------------------------------------------------------------- 15 16 !! * Modules used … … 26 27 27 28 !! * Routine accessibility 28 PUBLIC dyn_vor_enstrophy ! routine called by step.F90 29 PUBLIC dyn_vor_energy ! routine called by step.F90 30 PUBLIC dyn_vor_mixed ! routine called by step.F90 31 PUBLIC dyn_vor_ene_ens ! routine called by step.F90 32 PUBLIC dyn_vor_ctl ! routine called by step.F90 29 PUBLIC dyn_vor ! routine called by step.F90 33 30 34 31 !! * Shared module variables … … 38 35 LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme 39 36 37 !! * module variables 38 INTEGER :: & 39 nvor = 0 ! type of vorticity trend used 40 40 41 !! * Substitutions 41 42 # include "domzgr_substitute.h90" … … 43 44 !!---------------------------------------------------------------------- 44 45 !! OPA 9.0 , LOCEAN-IPSL (2005) 45 !! $Header$46 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt47 46 !!---------------------------------------------------------------------- 48 47 49 48 CONTAINS 50 49 51 SUBROUTINE dyn_vor_energy( kt ) 52 !!---------------------------------------------------------------------- 53 !! *** ROUTINE dyn_vor_energy *** 50 SUBROUTINE dyn_vor( kt ) 51 !!---------------------------------------------------------------------- 52 !! 53 !! ** Purpose : compute the lateral ocean tracer physics. 54 !! 55 !! ** Action : - Update (ua,va) with the now vorticity term trend 56 !! - save the trends in (utrd,vtrd) in 2 parts (relative 57 !! and planetary vorticity trends) ('key_trddyn') 58 !! 59 !! History : 60 !! 9.0 ! 05-11 (G. Madec) Original code 61 !!---------------------------------------------------------------------- 62 USE oce, ONLY : ztrdu => ta, & ! use ta as 3D workspace 63 ztrdv => sa ! use sa as 3D workspace 64 65 !! * Arguments 66 INTEGER, INTENT( in ) :: kt ! ocean time-step index 67 !!---------------------------------------------------------------------- 68 69 IF( kt == nit000 ) CALL vor_ctl ! initialisation & control of options 70 71 ! ! vorticity term including Coriolis 72 SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend 73 74 CASE ( -1 ) ! esopa: test all possibility with control print 75 CALL vor_ene( kt, 'TOT', ua, va ) 76 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 77 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 78 CALL vor_ens( kt, 'TOT', ua, va ) 79 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 80 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 81 CALL vor_mix( kt ) 82 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 83 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 84 CALL vor_een( kt, 'TOT', ua, va ) 85 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 86 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 87 88 CASE ( 0 ) ! energy conserving scheme 89 IF( l_trddyn ) THEN 90 ztrdu(:,:,:) = ua(:,:,:) 91 ztrdv(:,:,:) = va(:,:,:) 92 CALL vor_ene( kt, 'VOR', ua, va ) ! relative vorticity trend 93 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 94 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 95 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 96 ztrdu(:,:,:) = ua(:,:,:) 97 ztrdv(:,:,:) = va(:,:,:) 98 CALL vor_ene( kt, 'COR', ua, va ) ! planetary vorticity trend 99 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 102 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 103 ELSE 104 CALL vor_ene( kt, 'TOT', ua, va ) ! total vorticity 105 ENDIF 106 107 CASE ( 1 ) ! enstrophy conserving scheme 108 IF( l_trddyn ) THEN 109 ztrdu(:,:,:) = ua(:,:,:) 110 ztrdv(:,:,:) = va(:,:,:) 111 CALL vor_ens( kt, 'VOR', ua, va ) ! relative vorticity trend 112 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 113 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 114 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 115 ztrdu(:,:,:) = ua(:,:,:) 116 ztrdv(:,:,:) = va(:,:,:) 117 CALL vor_ens( kt, 'COR', ua, va ) ! planetary vorticity trend 118 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 119 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 120 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 121 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 122 ELSE 123 CALL vor_ens( kt, 'TOT', ua, va ) ! total vorticity 124 ENDIF 125 126 CASE ( 2 ) ! mixed ene-ens scheme 127 IF( l_trddyn ) THEN 128 ztrdu(:,:,:) = ua(:,:,:) 129 ztrdv(:,:,:) = va(:,:,:) 130 CALL vor_ens( kt, 'VOR', ua, va ) ! relative vorticity trend (ens) 131 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 134 ztrdu(:,:,:) = ua(:,:,:) 135 ztrdv(:,:,:) = va(:,:,:) 136 CALL vor_ene( kt, 'COR', ua, va ) ! planetary vorticity trend (ene) 137 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 138 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 139 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 140 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 141 ELSE 142 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) 143 ENDIF 144 145 CASE ( 3 ) ! energy and enstrophy conserving scheme 146 IF( l_trddyn ) THEN 147 ztrdu(:,:,:) = ua(:,:,:) 148 ztrdv(:,:,:) = va(:,:,:) 149 CALL vor_een( kt, 'VOR', ua, va ) ! relative vorticity trend 150 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 151 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 152 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 153 ztrdu(:,:,:) = ua(:,:,:) 154 ztrdv(:,:,:) = va(:,:,:) 155 CALL vor_een( kt, 'COR', ua, va ) ! planetary vorticity trend 156 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 157 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 158 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 159 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 160 ELSE 161 CALL vor_een( kt, 'TOT', ua, va ) ! total vorticity 162 ENDIF 163 164 END SELECT 165 166 ! ! print sum trends (used for debugging) 167 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 168 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 169 170 END SUBROUTINE dyn_vor 171 172 173 SUBROUTINE vor_ene( kt, cd_vor, pua, pva ) 174 !!---------------------------------------------------------------------- 175 !! *** ROUTINE vor_ene *** 54 176 !! 55 177 !! ** Purpose : Compute the now total vorticity trend and add it to … … 60 182 !! horizontal kinetic energy. 61 183 !! The trend of the vorticity term is given by: 62 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:184 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 63 185 !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] 64 186 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] … … 81 203 !! 9.0 ! 04-08 (C. Talandier) New trends organization 82 204 !!---------------------------------------------------------------------- 83 !! * Modules used84 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace85 ztdva => sa ! use sa as 3D workspace86 87 205 !! * Arguments 88 206 INTEGER, INTENT( in ) :: kt ! ocean time-step index 207 CHARACTER(len=3) , INTENT( in ) :: & 208 cd_vor ! define the vorticity considered 209 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 210 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 211 pua, pva ! ???? 89 212 90 213 !! * Local declarations 91 214 INTEGER :: ji, jj, jk ! dummy loop indices 92 215 REAL(wp) :: & 93 zfact2, zua, zva,& ! temporary scalars216 zfact2, & ! temporary scalars 94 217 zx1, zx2, zy1, zy2 ! " " 95 218 REAL(wp), DIMENSION(jpi,jpj) :: & 96 219 zwx, zwy, zwz ! temporary workspace 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &98 zcu, zcv ! " "99 220 !!---------------------------------------------------------------------- 100 221 101 222 IF( kt == nit000 ) THEN 102 223 IF(lwp) WRITE(numout,*) 103 IF(lwp) WRITE(numout,*) 'dyn _vor_energy: vorticity term: energy conserving scheme'104 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~'224 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 225 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 105 226 ENDIF 106 227 107 228 ! Local constant initialization 108 229 zfact2 = 0.5 * 0.5 109 zcu (:,:,:) = 0.e0 110 zcv (:,:,:) = 0.e0 111 112 ! Save ua and va trends 113 IF( l_trddyn ) THEN 114 ztdua(:,:,:) = ua(:,:,:) 115 ztdva(:,:,:) = va(:,:,:) 116 zcu(:,:,:) = 0.e0 117 zcv(:,:,:) = 0.e0 118 ENDIF 119 230 231 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 232 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 120 233 ! ! =============== 121 234 DO jk = 1, jpkm1 ! Horizontal slab … … 124 237 ! Potential vorticity and horizontal fluxes 125 238 ! ----------------------------------------- 126 IF( lk_sco ) THEN 127 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 239 SELECT CASE( cd_vor ) ! vorticity considered 240 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 241 zwz(:,:) = ff(:,:) 242 CASE ( 'VOR' ) ! relative vorticity 243 zwz(:,:) = rotn(:,:,jk) 244 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 245 zwz(:,:) = rotn(:,:,jk) + ff(:,:) 246 END SELECT 247 248 IF( ln_sco ) THEN 249 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 128 250 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 129 251 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 130 252 ELSE 131 zwz(:,:) = rotn(:,:,jk) + ff(:,:)132 253 zwx(:,:) = e2u(:,:) * un(:,:,jk) 133 254 zwy(:,:) = e1v(:,:) * vn(:,:,jk) … … 142 263 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 143 264 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 144 zua = zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 145 zva =-zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 146 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 147 va(ji,jj,jk) = va(ji,jj,jk) + zva 265 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 266 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 148 267 END DO 149 268 END DO … … 152 271 ! ! =============== 153 272 154 ! save the relative & planetary vorticity trends for diagnostic 155 ! momentum trends 156 IF( l_trddyn ) THEN 157 ! Compute the planetary vorticity term trend 158 ! ! =============== 159 DO jk = 1, jpkm1 ! Horizontal slab 160 ! ! =============== 161 DO jj = 2, jpjm1 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 164 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 165 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 166 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 167 # if defined key_s_coord 168 zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji ,jj-1) / fse3f(ji,jj-1,jk) * zy1 & 169 & + ff(ji ,jj ) / fse3f(ji,jj ,jk) * zy2 ) 170 zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj ) / fse3f(ji-1,jj,jk) * zx1 & 171 & + ff(ji ,jj ) / fse3f(ji ,jj,jk) * zx2 ) 172 # else 173 zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 174 zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 175 # endif 176 END DO 177 END DO 178 ! ! =============== 179 END DO ! End of slab 180 ! ! =============== 181 182 ! Compute the relative vorticity term trend 183 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 184 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 185 186 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 187 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 188 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 189 ENDIF 190 191 IF(ln_ctl) THEN ! print sum trends (used for debugging) 192 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 193 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 194 ENDIF 195 196 END SUBROUTINE dyn_vor_energy 197 198 199 SUBROUTINE dyn_vor_mixed( kt ) 200 !!---------------------------------------------------------------------- 201 !! *** ROUTINE dyn_vor_mixed *** 273 END SUBROUTINE vor_ene 274 275 276 SUBROUTINE vor_mix( kt ) 277 !!---------------------------------------------------------------------- 278 !! *** ROUTINE vor_mix *** 202 279 !! 203 280 !! ** Purpose : Compute the now total vorticity trend and add it to … … 209 286 !! ticity term and the horizontal kinetic energy for (f x uh), the 210 287 !! coriolis term. the now trend of the vorticity term is given by: 211 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:288 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 212 289 !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] 213 290 !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ] … … 234 311 !! 9.0 ! 04-08 (C. Talandier) New trends organization 235 312 !!---------------------------------------------------------------------- 236 !! * Modules used237 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace238 ztdva => sa ! use sa as 3D workspace239 313 !! * Arguments 240 314 INTEGER, INTENT( in ) :: kt ! ocean timestep index … … 247 321 REAL(wp), DIMENSION(jpi,jpj) :: & 248 322 zwx, zwy, zwz, zww ! temporary workspace 249 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &250 zcu, zcv ! " "251 323 !!---------------------------------------------------------------------- 252 324 253 325 IF( kt == nit000 ) THEN 254 326 IF(lwp) WRITE(numout,*) 255 IF(lwp) WRITE(numout,*) 'dyn _vor_mixed: vorticity term: mixed energy/enstrophy conserving scheme'256 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~'327 IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' 328 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 257 329 ENDIF 258 330 … … 261 333 zfact2 = 0.5 * 0.5 262 334 263 ! Save ua and va trends 264 IF( l_trddyn ) THEN 265 ztdua(:,:,:) = ua(:,:,:) 266 ztdva(:,:,:) = va(:,:,:) 267 zcu(:,:,:) = 0.e0 268 zcv(:,:,:) = 0.e0 269 ENDIF 270 335 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 336 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 271 337 ! ! =============== 272 338 DO jk = 1, jpkm1 ! Horizontal slab … … 275 341 ! Relative and planetary potential vorticity and horizontal fluxes 276 342 ! ---------------------------------------------------------------- 277 IF( l k_sco ) THEN343 IF( ln_sco ) THEN 278 344 zwz(:,:) = ff (:,:) / fse3f(:,:,jk) 279 345 zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) … … 310 376 ! ! =============== 311 377 312 ! save the relative & planetary vorticity trends for diagnostic 313 ! momentum trends 314 IF( l_trddyn ) THEN 315 ! Compute the planetary vorticity term trend 316 ! ! =============== 317 DO jk = 1, jpkm1 ! Horizontal slab 318 ! ! =============== 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 322 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 323 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 324 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 325 ! energy conserving formulation for planetary vorticity term 326 zcu(ji,jj,jk) = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 327 zcv(ji,jj,jk) =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 328 END DO 329 END DO 330 ! ! =============== 331 END DO ! End of slab 332 ! ! =============== 333 334 ! Compute the relative vorticity term trend 335 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 336 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 337 338 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 339 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 340 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 341 ENDIF 342 343 IF(ln_ctl) THEN ! print sum trends (used for debugging) 344 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 345 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 346 ENDIF 347 348 END SUBROUTINE dyn_vor_mixed 349 350 351 SUBROUTINE dyn_vor_enstrophy( kt ) 352 !!---------------------------------------------------------------------- 353 !! *** ROUTINE dyn_vor_enstrophy *** 378 END SUBROUTINE vor_mix 379 380 381 SUBROUTINE vor_ens( kt, cd_vor, pua, pva ) 382 !!---------------------------------------------------------------------- 383 !! *** ROUTINE vor_ens *** 354 384 !! 355 385 !! ** Purpose : Compute the now total vorticity trend and add it to … … 360 390 !! potential enstrophy of a horizontally non-divergent flow. the 361 391 !! trend of the vorticity term is given by: 362 !! * s-coordinate (l k_sco=T), the e3. are inside the derivative:392 !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: 363 393 !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 364 394 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] … … 381 411 !! 9.0 ! 04-08 (C. Talandier) New trends organization 382 412 !!---------------------------------------------------------------------- 383 !! * modules used384 USE oce, ONLY: zwx => ta, & ! use ta as 3D workspace385 zwy => sa ! use sa as 3D workspace386 413 !! * Arguments 387 414 INTEGER, INTENT( in ) :: kt ! ocean timestep 415 CHARACTER(len=3) , INTENT( in ) :: & 416 cd_vor ! define the vorticity considered 417 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 418 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 419 pua, pva ! ???? 388 420 389 421 !! * Local declarations 390 422 INTEGER :: ji, jj, jk ! dummy loop indices 391 423 REAL(wp) :: & 392 zfact1, zua, zva, zuav, zvau ! temporary scalars 393 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 394 zcu, zcv, zwz, & ! temporary workspace 395 ztdua, ztdva ! temporary workspace 424 zfact1, zuav, zvau ! temporary scalars 425 REAL(wp), DIMENSION(jpi,jpj) :: & 426 zwx, zwy, zwz ! temporary workspace 396 427 !!---------------------------------------------------------------------- 397 428 398 429 IF( kt == nit000 ) THEN 399 430 IF(lwp) WRITE(numout,*) 400 IF(lwp) WRITE(numout,*) 'dyn _vor_enstrophy: vorticity term: enstrophy conserving scheme'401 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~~~~'431 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 432 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 402 433 ENDIF 403 434 … … 405 436 zfact1 = 0.5 * 0.25 406 437 407 ! Save ua and va trends 408 IF( l_trddyn ) THEN 409 ztdua(:,:,:) = ua(:,:,:) 410 ztdva(:,:,:) = va(:,:,:) 411 zcu(:,:,:) = 0.e0 412 zcv(:,:,:) = 0.e0 413 ENDIF 414 438 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 439 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 415 440 ! ! =============== 416 441 DO jk = 1, jpkm1 ! Horizontal slab … … 419 444 ! Potential vorticity and horizontal fluxes 420 445 ! ----------------------------------------- 421 IF( lk_sco ) THEN 446 SELECT CASE( cd_vor ) ! vorticity considered 447 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 448 zwz(:,:) = ff(:,:) 449 CASE ( 'VOR' ) ! relative vorticity 450 zwz(:,:) = rotn(:,:,jk) 451 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 452 zwz(:,:) = rotn(:,:,jk) + ff(:,:) 453 END SELECT 454 455 IF( ln_sco ) THEN 422 456 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 423 457 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 424 zwz(ji,jj ,jk) = ( rotn(ji,jj,jk) + ff(ji,jj)) / fse3f(ji,jj,jk)425 zwx(ji,jj ,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)426 zwy(ji,jj ,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)458 zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 459 zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 460 zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 427 461 END DO 428 462 END DO … … 430 464 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 431 465 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 432 zwz(ji,jj,jk) = rotn(ji,jj,jk) + ff(ji,jj) 433 zwx(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 434 zwy(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 466 zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 467 zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 435 468 END DO 436 469 END DO … … 442 475 DO jj = 2, jpjm1 443 476 DO ji = fs_2, fs_jpim1 ! vector opt. 444 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1,jk) + zwy(ji+1,jj-1,jk) & 445 + zwy(ji ,jj ,jk) + zwy(ji+1,jj ,jk) ) 446 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ,jk) + zwx(ji-1,jj+1,jk) & 447 + zwx(ji ,jj ,jk) + zwx(ji ,jj+1,jk) ) 448 449 zua = zuav * ( zwz(ji ,jj-1,jk) + zwz(ji,jj,jk) ) 450 zva = zvau * ( zwz(ji-1,jj ,jk) + zwz(ji,jj,jk) ) 451 452 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 453 va(ji,jj,jk) = va(ji,jj,jk) + zva 477 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 478 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 479 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 480 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 481 482 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 483 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 454 484 END DO 455 485 END DO … … 458 488 ! ! =============== 459 489 460 461 ! save the relative & planetary vorticity trends for diagnostic 462 ! momentum trends 463 IF( l_trddyn ) THEN 464 ! Compute the planetary vorticity term trend 465 ! ! =============== 466 DO jk = 1, jpkm1 ! Horizontal slab 467 ! ! =============== 468 DO jj = 2, jpjm1 469 DO ji = fs_2, fs_jpim1 ! vector opt. 470 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1,jk) + zwy(ji+1,jj-1,jk) & 471 & + zwy(ji ,jj ,jk) + zwy(ji+1,jj ,jk) ) 472 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ,jk) + zwx(ji-1,jj+1,jk) & 473 & + zwx(ji ,jj ,jk) + zwx(ji ,jj+1,jk) ) 474 # if defined key_s_coord 475 zcu(ji,jj,jk) = zuav * ( ff(ji ,jj-1) / fse3f(ji ,jj-1,jk) & 476 & + ff(ji ,jj ) / fse3f(ji ,jj ,jk) ) 477 zcv(ji,jj,jk) = zvau * ( ff(ji-1,jj ) / fse3f(ji-1,jj ,jk) & 478 & + ff(ji ,jj ) / fse3f(ji ,jj ,jk) ) 479 # else 480 zcu(ji,jj,jk) = zuav * ( ff(ji ,jj-1) + ff(ji,jj) ) 481 zcv(ji,jj,jk) = zvau * ( ff(ji-1,jj ) + ff(ji,jj) ) 482 # endif 483 END DO 484 END DO 485 ! ! =============== 486 END DO ! End of slab 487 ! ! =============== 488 489 ! Compute the relative vorticity term trend 490 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 491 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 492 493 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 494 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 495 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 496 ENDIF 497 498 IF(ln_ctl) THEN ! print sum trends (used for debugging) 499 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 500 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 501 ENDIF 502 503 END SUBROUTINE dyn_vor_enstrophy 504 505 506 SUBROUTINE dyn_vor_ene_ens( kt ) 507 !!---------------------------------------------------------------------- 508 !! *** ROUTINE dyn_vor_ene_ens *** 490 END SUBROUTINE vor_ens 491 492 493 SUBROUTINE vor_een( kt, cd_vor, pua, pva ) 494 !!---------------------------------------------------------------------- 495 !! *** ROUTINE vor_een *** 509 496 !! 510 497 !! ** Purpose : Compute the now total vorticity trend and add it to … … 516 503 !! when horizontal divergence is zero. 517 504 !! The trend of the vorticity term is given by: 518 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:505 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 519 506 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 520 507 !! Add this trend to the general momentum trend (ua,va): … … 534 521 !! 9.0 ! 04-08 (C. Talandier) New trends organization 535 522 !!---------------------------------------------------------------------- 536 !! * Modules used537 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace538 ztdva => sa ! use sa as 3D workspace539 540 523 !! * Arguments 541 524 INTEGER, INTENT( in ) :: kt ! ocean time-step index 525 CHARACTER(len=3) , INTENT( in ) :: & 526 cd_vor ! define the vorticity considered 527 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 528 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 529 pua, pva ! ???? 542 530 543 531 !! * Local declarations … … 547 535 REAL(wp), DIMENSION(jpi,jpj) :: & 548 536 zwx, zwy, zwz, & ! temporary workspace 549 ztnw, ztne, ztsw, ztse, & ! " " 550 zcor ! potential planetary vorticity (f/e3) 551 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 552 zcu, zcv ! temporary workspace 537 ztnw, ztne, ztsw, ztse ! " " 553 538 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: & 554 539 ze3f … … 557 542 IF( kt == nit000 ) THEN 558 543 IF(lwp) WRITE(numout,*) 559 IF(lwp) WRITE(numout,*) 'dyn _vor_ene_ens: vorticity term: energy and enstrophy conserving scheme'560 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~~'544 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 545 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 561 546 562 547 DO jk = 1, jpk … … 565 550 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 566 551 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25 567 !!! ze3f(ji,jj,jk) = MAX( ze3f(ji,jj,jk) , 1.e-20)568 552 IF( ze3f(ji,jj,jk) /= 0.e0 ) ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk) 569 553 END DO … … 576 560 zfac12 = 1.e0 / 12.e0 577 561 578 ! Save ua and va trends579 IF( l_trddyn ) THEN580 ztdua(:,:,:) = ua(:,:,:)581 ztdva(:,:,:) = va(:,:,:)582 zcu(:,:,:) = 0.e0583 zcv(:,:,:) = 0.e0584 ENDIF585 562 563 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 564 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 586 565 ! ! =============== 587 566 DO jk = 1, jpkm1 ! Horizontal slab … … 590 569 ! Potential vorticity and horizontal fluxes 591 570 ! ----------------------------------------- 592 !!!bug zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 593 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 571 SELECT CASE( cd_vor ) ! vorticity considered 572 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 573 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 574 CASE ( 'VOR' ) ! relative vorticity 575 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 576 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 577 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 578 END SELECT 579 594 580 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 595 581 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 596 zcor(:,:) = ff(:,:) * ze3f(:,:,jk)597 582 598 583 ! Compute and add the vorticity term trend … … 621 606 zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 622 607 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 623 ua(ji,jj,jk) =ua(ji,jj,jk) + zua624 va(ji,jj,jk) =va(ji,jj,jk) + zva608 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 609 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 625 610 END DO 626 611 END DO … … 629 614 ! ! =============== 630 615 631 ! save the relative & planetary vorticity trends for diagnostic 632 ! momentum trends 633 IF( l_trddyn ) THEN 634 ! Compute the planetary vorticity term trend 635 ! ! =============== 636 DO jk = 1, jpkm1 ! Horizontal slab 637 ! ! =============== 638 DO jj = 2, jpjm1 639 DO ji = fs_2, fs_jpim1 ! vector opt. 640 zcu(ji,jj,jk) = + zfac12 / e1u(ji,jj) * ( zcor(ji,jj ) * zwy(ji ,jj ) + zcor(ji+1,jj) * zwy(ji+1,jj ) & 641 & + zcor(ji,jj ) * zwy(ji ,jj-1) + zcor(ji+1,jj) * zwy(ji+1,jj-1) ) 642 zcv(ji,jj,jk) = - zfac12 / e2v(ji,jj) * ( zcor(ji,jj+1) * zwx(ji-1,jj+1) + zcor(ji,jj+1) * zwx(ji ,jj+1) & 643 & + zcor(ji,jj ) * zwx(ji-1,jj ) + zcor(ji,jj ) * zwx(ji ,jj ) ) 644 END DO 645 END DO 646 ! ! =============== 647 END DO ! End of slab 648 ! ! =============== 649 650 ! Compute the relative vorticity term trend 651 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 652 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 653 654 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 655 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 656 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 657 ENDIF 658 659 IF(ln_ctl) THEN ! print sum trends (used for debugging) 660 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 661 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 662 ENDIF 663 664 END SUBROUTINE dyn_vor_ene_ens 665 666 667 SUBROUTINE dyn_vor_ctl 616 END SUBROUTINE vor_een 617 618 619 SUBROUTINE vor_ctl 668 620 !!--------------------------------------------------------------------- 669 !! *** ROUTINE dyn_vor_ctl ***621 !! *** ROUTINE vor_ctl *** 670 622 !! 671 623 !! ** Purpose : Control the consistency between cpp options for … … 691 643 IF(lwp) THEN 692 644 WRITE(numout,*) 693 WRITE(numout,*) 'dyn _vor_ctl : vorticity term : read namelist and control the consistency'645 WRITE(numout,*) 'dyn:vor_ctl : vorticity term : read namelist and control the consistency' 694 646 WRITE(numout,*) '~~~~~~~~~~~' 695 647 WRITE(numout,*) ' Namelist nam_dynvor : oice of the vorticity term scheme' 648 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 696 649 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 697 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene698 650 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 699 651 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een … … 701 653 702 654 ioptio = 0 703 655 IF( ln_dynvor_ene ) THEN 656 nvor = 0 657 IF(lwp) WRITE(numout,*) 658 IF(lwp) WRITE(numout,*) ' vorticity term : energy conserving scheme' 659 ioptio = ioptio + 1 660 ENDIF 704 661 IF( ln_dynvor_ens ) THEN 662 nvor = 1 705 663 IF(lwp) WRITE(numout,*) 706 664 IF(lwp) WRITE(numout,*) ' vorticity term : enstrophy conserving scheme' 707 665 ioptio = ioptio + 1 708 666 ENDIF 709 IF( ln_dynvor_ene ) THEN710 IF(lwp) WRITE(numout,*)711 IF(lwp) WRITE(numout,*) ' vorticity term : energy conserving scheme'712 ioptio = ioptio + 1713 ENDIF714 667 IF( ln_dynvor_mix ) THEN 668 nvor = 2 715 669 IF(lwp) WRITE(numout,*) 716 670 IF(lwp) WRITE(numout,*) ' vorticity term : mixed enstrophy/energy conserving scheme' … … 718 672 ENDIF 719 673 IF( ln_dynvor_een ) THEN 674 nvor = 3 720 675 IF(lwp) WRITE(numout,*) 721 676 IF(lwp) WRITE(numout,*) ' vorticity term : energy and enstrophy conserving scheme' 722 677 ioptio = ioptio + 1 723 678 ENDIF 724 IF ( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 725 WRITE(numout,cform_err) 726 IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 727 nstop = nstop + 1 728 ENDIF 729 730 END SUBROUTINE dyn_vor_ctl 679 IF( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 680 WRITE(numout,cform_err) 681 IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 682 nstop = nstop + 1 683 ENDIF 684 IF( lk_esopa ) THEN 685 nvor = -1 686 IF(lwp ) WRITE(numout,*) ' esopa test: use all lateral physics options' 687 ENDIF 688 IF(lwp) WRITE(numout,*) ' choice of vor_... nvor = ', nvor 689 690 END SUBROUTINE vor_ctl 731 691 732 692 !!============================================================================== -
trunk/NEMO/OPA_SRC/DYN/dynzad.F90
r258 r455 34 34 CONTAINS 35 35 36 #if defined key_ autotasking37 !!---------------------------------------------------------------------- 38 !! 'key_ autotasking'j-k-i loops (j-slab)36 #if defined key_mpp_omp 37 !!---------------------------------------------------------------------- 38 !! 'key_mpp_omp' OpenMP / NEC autotasking: j-k-i loops (j-slab) 39 39 !!---------------------------------------------------------------------- 40 40 … … 46 46 !! add it to the general trend of momentum equation. 47 47 !! 48 !! ** Method : Use j-slab (j-k-i loops) for auto-tasking48 !! ** Method : Use j-slab (j-k-i loops) for OpenMP / NEC autotasking 49 49 !! The now vertical advection of momentum is given by: 50 50 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] -
trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r258 r455 16 16 USE in_out_manager ! I/O manager 17 17 USE taumod ! surface ocean stress 18 USE trdmod ! ocean dynamics trends19 USE trdmod_oce ! ocean variables trends20 18 USE prtctl ! Print control 21 19 … … 68 66 !! * Local declarations 69 67 INTEGER :: & 70 ji, jj, jk, jl, & ! dummy loop indices 71 ikbu, ikbum1 , ikbv, ikbvm1 ! temporary integers 68 ji, jj, jk, jl ! dummy loop indices 72 69 REAL(wp) :: & 73 70 zrau0r, zlavmr, z2dt, zua, zva ! temporary scalars 74 71 REAL(wp), DIMENSION(jpi,jpk) :: & 75 72 zwx, zwy, zwz, zww ! temporary workspace arrays 76 REAL(wp), DIMENSION(jpi,jpj) :: &77 ztsx, ztsy, ztbx, ztby ! temporary workspace arrays78 73 !!---------------------------------------------------------------------- 79 74 … … 89 84 zlavmr = 1. / float( n_zdfexp ) ! inverse of the number of sub time step 90 85 z2dt = 2. * rdt ! Leap-frog environnement 91 ztsx(:,:) = 0.e092 ztsy(:,:) = 0.e093 ztbx(:,:) = 0.e094 ztby(:,:) = 0.e095 96 ! Save ua and va trends97 IF( l_trddyn ) THEN98 ztdua(:,:,:) = ua(:,:,:)99 ztdva(:,:,:) = va(:,:,:)100 ENDIF101 86 102 87 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! Euler time stepping when starting from rest … … 150 135 ! ! =============== 151 136 152 ! save the vertical diffusion trends for diagnostic153 ! momentum trends154 IF( l_trddyn ) THEN155 ! save the total vertical momentum diffusive trend156 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)157 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)158 159 ! subtract and save surface and momentum fluxes160 ! ! ===============161 DO jj = 2, jpjm1 ! Horizontal slab162 ! ! ===============163 DO ji = 2, jpim1164 ! save the surface momentum fluxes165 ztsx(ji,jj) = zwy(ji,1) / fse3u(ji,jj,1)166 ztsy(ji,jj) = zww(ji,1) / fse3v(ji,jj,1)167 ! save bottom friction momentum fluxes168 ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )169 ikbum1 = MAX( ikbu-1, 1 )170 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )171 ikbvm1 = MAX( ikbv-1, 1 )172 ztbx(ji,jj) = avmu(ji,jj,ikbu) * zwx(ji,ikbum1) &173 / ( fse3u(ji,jj,ikbum1) * fse3uw(ji,jj,ikbu) )174 ztby(ji,jj) = avmv(ji,jj,ikbv) * zwz(ji,ikbvm1) &175 / ( fse3v(ji,jj,ikbvm1) * fse3vw(ji,jj,ikbv) )176 ! subtract surface forcing and bottom friction trend from vertical177 ! diffusive momentum trend178 ztdua(ji,jj,1 ) = ztdua(ji,jj,1 ) - ztsx(ji,jj)179 ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj)180 ztdva(ji,jj,1 ) = ztdva(ji,jj,1 ) - ztsy(ji,jj)181 ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)182 END DO183 ! ! ===============184 END DO ! End of slab185 ! ! ===============186 187 CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt)188 ztdua(:,:,:) = 0.e0189 ztdva(:,:,:) = 0.e0190 ztdua(:,:,1) = ztsx(:,:)191 ztdva(:,:,1) = ztsy(:,:)192 CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt)193 ztdua(:,:,:) = 0.e0194 ztdva(:,:,:) = 0.e0195 ztdua(:,:,1) = ztbx(:,:)196 ztdva(:,:,1) = ztby(:,:)197 CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt)198 ENDIF199 200 IF(ln_ctl) THEN ! print sum trends (used for debugging)201 CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, &202 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')203 ENDIF204 205 137 END SUBROUTINE dyn_zdf_exp 206 138 -
trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r258 r455 20 20 USE in_out_manager ! I/O manager 21 21 USE taumod ! surface ocean stress 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 22 USE prtctl ! Print control 25 23 … … 60 58 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 61 59 !! mixing trend. 62 !! - Save the trends in (ztdua,ztdva) ('l_trddyn')63 60 !! 64 61 !! History : … … 76 73 77 74 !! * Local declarations 78 INTEGER :: & 79 ji, jj, jk, & ! dummy loop indices 80 ikbu, ikbum1, ikbv, ikbvm1 ! temporary integers 75 INTEGER :: ji, jj, jk ! dummy loop indices 81 76 REAL(wp) :: & 82 77 zrau0r, z2dt, & ! temporary scalars 83 78 z2dtf, zcoef, zzws, zrhs ! " " 84 REAL(wp), DIMENSION(jpi,jpj) :: &85 ztsx, ztsy, ztbx, ztby ! temporary workspace arrays86 79 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 87 zwi , ztdua, ztdva! temporary workspace arrays80 zwi ! temporary workspace arrays 88 81 !!---------------------------------------------------------------------- 89 82 … … 98 91 zrau0r = 1. / rau0 ! inverse of the reference density 99 92 z2dt = 2. * rdt ! Leap-frog environnement 100 ztsx(:,:) = 0.e0 101 ztsy(:,:) = 0.e0 102 ztbx(:,:) = 0.e0 103 ztby(:,:) = 0.e0 93 104 94 ! Euler time stepping when starting from rest 105 95 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 106 107 ! Save previous ua and va trends108 IF( l_trddyn ) THEN109 ztdua(:,:,:) = ua(:,:,:)110 ztdva(:,:,:) = va(:,:,:)111 ENDIF112 96 113 97 ! 1. Vertical diffusion on u … … 194 178 END DO 195 179 END DO 196 197 IF( l_trddyn ) THEN198 ! diagnose surface and bottom momentum fluxes199 DO jj = 2, jpjm1200 DO ji = fs_2, fs_jpim1 ! vector opt.201 ! save the surface forcing momentum fluxes202 ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )203 ! save bottom friction momentum fluxes204 ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )205 ikbum1 = MAX( ikbu-1, 1 )206 ztbx(ji,jj) = - avmu(ji,jj,ikbu) * ua(ji,jj,ikbum1) &207 / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )208 ! subtract surface forcing and bottom friction trend from vertical209 ! diffusive momentum trend210 ztdua(ji,jj,1 ) = ztdua(ji,jj,1 ) - ztsx(ji,jj)211 ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj)212 END DO213 END DO214 ENDIF215 180 216 181 ! Normalization to obtain the general momentum trend ua … … 309 274 END DO 310 275 311 IF( l_trddyn ) THEN 312 ! diagnose surface and bottom momentum fluxes 313 DO jj = 2, jpjm1 314 DO ji = fs_2, fs_jpim1 ! vector opt. 315 ! save the surface forcing momentum fluxes 316 ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 317 ! save bottom friction momentum fluxes 318 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 319 ikbvm1 = MAX( ikbv-1, 1 ) 320 ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1) & 321 / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 322 ! subtract surface forcing and bottom friction trend from vertical 323 ! diffusive momentum trend 324 ztdva(ji,jj,1 ) = ztdva(ji,jj,1 ) - ztsy(ji,jj) 325 ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj) 326 END DO 327 END DO 328 ENDIF 276 ! flux de surface doit etre calcule dans trdmod et boootom stress 277 ! deduit par integration verticale dans trmod pour jpdtdzdf 278 !RB IF( l_trddyn ) THEN 279 ! ! diagnose surface and bottom momentum fluxes 280 ! DO jj = 2, jpjm1 281 ! DO ji = fs_2, fs_jpim1 ! vector opt. 282 ! ! save the surface forcing momentum fluxes 283 ! ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 284 ! ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) 285 ! ! save bottom friction momentum fluxes 286 ! ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 287 ! ikbvm1 = MAX( ikbv-1, 1 ) 288 ! ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1) & 289 ! / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 290 ! ! subtract surface forcing and bottom friction trend from vertical 291 ! ! diffusive momentum trend 292 ! ztdva(ji,jj,1 ) = ztdva(ji,jj,1 ) - ztsy(ji,jj) 293 ! ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj) 294 ! END DO 295 ! END DO 296 ! ENDIF 329 297 330 298 ! Normalization to obtain the general momentum trend va … … 337 305 END DO 338 306 339 ! save the vertical diffusion trends for diagnostic340 ! momentum trends341 IF( l_trddyn ) THEN342 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)343 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)344 345 CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt)346 ztdua(:,:,:) = 0.e0347 ztdva(:,:,:) = 0.e0348 ztdua(:,:,1) = ztsx(:,:)349 ztdva(:,:,1) = ztsy(:,:)350 CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt)351 ztdua(:,:,:) = 0.e0352 ztdva(:,:,:) = 0.e0353 ztdua(:,:,1) = ztbx(:,:)354 ztdva(:,:,1) = ztby(:,:)355 CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt)356 ENDIF357 358 IF(ln_ctl) THEN ! print sum trends (used for debugging)359 CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, &360 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')361 ENDIF362 363 307 END SUBROUTINE dyn_zdf_imp 364 308 -
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r258 r455 26 26 CONTAINS 27 27 28 #if defined key_ autotasking28 #if defined key_mpp_omp 29 29 !!---------------------------------------------------------------------- 30 !! 'key_ autotasking'j-k-i loop (j-slab)30 !! 'key_mpp_omp' j-k-i loop (j-slab) 31 31 !!---------------------------------------------------------------------- 32 32 … … 64 64 IF(lwp) WRITE(numout,*) 65 65 IF(lwp) WRITE(numout,*) 'wzv : vertical velocity from continuity eq.' 66 IF(lwp) WRITE(numout,*) '~~~~~~~ auto-tasking case : j-k-i loop'66 IF(lwp) WRITE(numout,*) '~~~~~~~ j-k-i loops' 67 67 68 68 ! bottom boundary condition: w=0 (set once for all)
Note: See TracChangeset
for help on using the changeset viewer.