Changeset 2148
- Timestamp:
- 2010-10-04T15:53:42+02:00 (14 years ago)
- Location:
- branches/DEV_r2106_LOCEAN2010/NEMO
- Files:
-
- 1 added
- 35 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/dom_oce.F90
r2006 r2148 145 145 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3vw_1 !: analytical vertical scale factors at VW-- 146 146 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3w_1 , e3uw_1 !: W--UW points (m) 147 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3t_b !: before - - - - T points (m) 148 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3u_b , e3v_b !: - - - - - U--V points (m) 147 149 #else 148 150 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/domvvl.F90
r2004 r2148 51 51 !!---------------------------------------------------------------------- 52 52 INTEGER :: ji, jj, jk 53 REAL(wp) :: zcoefu, zcoefv, zcoeff 53 REAL(wp) :: zcoefu , zcoefv , zcoeff ! temporary scalars 54 REAL(wp) :: zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1 ! - - 55 REAL(wp), DIMENSION(jpi,jpj) :: zs_t, zs_u_1, zs_v_1 ! - 2D workspace 54 56 !!---------------------------------------------------------------------- 55 57 … … 62 64 IF( lk_zco ) CALL ctl_stop( 'dom_vvl : key_zco is incompatible with variable volume option key_vvl') 63 65 64 IF( ln_zco) THEN 65 DO jk = 1, jpk 66 gdept(:,:,jk) = gdept_0(jk) 67 gdepw(:,:,jk) = gdepw_0(jk) 68 gdep3w(:,:,jk) = gdepw_0(jk) 69 e3t (:,:,jk) = e3t_0(jk) 70 e3u (:,:,jk) = e3t_0(jk) 71 e3v (:,:,jk) = e3t_0(jk) 72 e3f (:,:,jk) = e3t_0(jk) 73 e3w (:,:,jk) = e3w_0(jk) 74 e3uw(:,:,jk) = e3w_0(jk) 75 e3vw(:,:,jk) = e3w_0(jk) 76 END DO 77 ELSE 78 fsdept(:,:,:) = gdept (:,:,:) 79 fsdepw(:,:,:) = gdepw (:,:,:) 80 fsde3w(:,:,:) = gdep3w(:,:,:) 81 fse3t (:,:,:) = e3t (:,:,:) 82 fse3u (:,:,:) = e3u (:,:,:) 83 fse3v (:,:,:) = e3v (:,:,:) 84 fse3f (:,:,:) = e3f (:,:,:) 85 fse3w (:,:,:) = e3w (:,:,:) 86 fse3uw(:,:,:) = e3uw (:,:,:) 87 fse3vw(:,:,:) = e3vw (:,:,:) 88 ENDIF 66 fsdept(:,:,:) = gdept (:,:,:) 67 fsdepw(:,:,:) = gdepw (:,:,:) 68 fsde3w(:,:,:) = gdep3w(:,:,:) 69 fse3t (:,:,:) = e3t (:,:,:) 70 fse3u (:,:,:) = e3u (:,:,:) 71 fse3v (:,:,:) = e3v (:,:,:) 72 fse3f (:,:,:) = e3f (:,:,:) 73 fse3w (:,:,:) = e3w (:,:,:) 74 fse3uw(:,:,:) = e3uw (:,:,:) 75 fse3vw(:,:,:) = e3vw (:,:,:) 89 76 90 77 ! !== mu computation ==! … … 130 117 hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 131 118 END DO 119 120 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 121 ! for ssh and scale factors 122 zs_t (:,:) = e1t(:,:) * e2t(:,:) 123 zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 124 zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 132 125 133 126 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 134 127 DO ji = 1, jpim1 ! NO vector opt. 135 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 136 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 137 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 138 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 139 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 140 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 141 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 142 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & 143 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 144 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 145 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 146 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 147 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 148 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 149 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 128 zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj) 129 zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj) 130 zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) ) 131 ! before fields 132 zv_t_ij = zs_t(ji ,jj ) * sshb(ji ,jj ) 133 zv_t_ip1j = zs_t(ji+1,jj ) * sshb(ji+1,jj ) 134 zv_t_ijp1 = zs_t(ji ,jj+1) * sshb(ji ,jj+1) 135 sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 136 sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 137 ! now fields 138 zv_t_ij = zs_t(ji ,jj ) * sshn(ji ,jj ) 139 zv_t_ip1j = zs_t(ji+1,jj ) * sshn(ji+1,jj ) 140 zv_t_ijp1 = zs_t(ji ,jj+1) * sshn(ji ,jj+1) 141 zv_t_ip1jp1 = zs_t(ji ,jj+1) * sshn(ji ,jj+1) 142 sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 143 sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 144 sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 ) 150 145 END DO 151 146 END DO 152 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) ! lateral boundary conditions 153 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 154 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 147 CALL lbc_lnk( sshu_n, 'U', 1. ) ; CALL lbc_lnk( sshu_b, 'U', 1. ) ! lateral boundary conditions 148 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 149 CALL lbc_lnk( sshf_n, 'F', 1. ) 150 151 ! initialise before scale factors at (u/v)-points 152 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 153 DO jk = 1, jpkm1 154 DO jj = 1, jpjm1 155 DO ji = 1, jpim1 156 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 157 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 158 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 159 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 160 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 161 END DO 162 END DO 163 END DO 164 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 165 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 166 ! Add initial scale factor to scale factor anomaly 167 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 168 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 155 169 ! 156 DO jk = 1, jpkm1157 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays158 fsdepw(:,:,jk) = fsdepw_n(:,:,jk)159 fsde3w(:,:,jk) = fsde3w_n(:,:,jk)160 !161 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays162 fse3u (:,:,jk) = fse3u_n (:,:,jk)163 fse3v (:,:,jk) = fse3v_n (:,:,jk)164 fse3f (:,:,jk) = fse3f_n (:,:,jk)165 fse3w (:,:,jk) = fse3w_n (:,:,jk)166 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)167 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)168 END DO169 170 171 172 170 END SUBROUTINE dom_vvl 173 171 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
r1565 r2148 46 46 # define fse3vw(i,j,k) e3vw_1(i,j,k) 47 47 48 # define fsdept_b(i,j,k) (fsdept_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 49 # define fsdepw_b(i,j,k) (fsdepw_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 50 # define fsde3w_b(i,j,k) (fsde3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))-sshb(i,j)) 51 # define fse3t_b(i,j,k) (fse3t_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 52 # define fse3u_b(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 53 # define fse3v_b(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 54 # define fse3f_b(i,j,k) (fse3f_0(i,j,k)*(1.+sshf_b(i,j)*muf(i,j,k))) 55 # define fse3w_b(i,j,k) (fse3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 48 # define fse3t_b(i,j,k) e3t_b(i,j,k) 49 # define fse3u_b(i,j,k) e3u_b(i,j,k) 50 # define fse3v_b(i,j,k) e3v_b(i,j,k) 56 51 # define fse3uw_b(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 57 52 # define fse3vw_b(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) … … 70 65 # define fse3t_m(i,j,k) (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 71 66 72 # define fsdept_a(i,j,k) (fsdept_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))73 # define fsdepw_a(i,j,k) (fsdepw_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))74 # define fsde3w_a(i,j,k) (fsde3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))-ssha(i,j))75 67 # define fse3t_a(i,j,k) (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 76 68 # define fse3u_a(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 77 69 # define fse3v_a(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 78 # define fse3f_a(i,j,k) (fse3f_0(i,j,k)*(1.+sshf_a(i,j)*muf(i,j,k)))79 # define fse3w_a(i,j,k) (fse3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))80 # define fse3uw_a(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k)))81 # define fse3vw_a(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k)))82 70 83 71 #else … … 94 82 # define fse3vw(i,j,k) fse3vw_0(i,j,k) 95 83 96 # define fsdept_b(i,j,k) fsdept_0(i,j,k)97 # define fsdepw_b(i,j,k) fsdepw_0(i,j,k)98 # define fsde3w_b(i,j,k) fsde3w_0(i,j,k)99 84 # define fse3t_b(i,j,k) fse3t_0(i,j,k) 100 85 # define fse3u_b(i,j,k) fse3u_0(i,j,k) 101 86 # define fse3v_b(i,j,k) fse3v_0(i,j,k) 102 # define fse3f_b(i,j,k) fse3f_0(i,j,k)103 # define fse3w_b(i,j,k) fse3w_0(i,j,k)104 87 # define fse3uw_b(i,j,k) fse3uw_0(i,j,k) 105 88 # define fse3vw_b(i,j,k) fse3vw_0(i,j,k) … … 118 101 # define fse3t_m(i,j,k) fse3t_0(i,j,k) 119 102 120 # define fsdept_a(i,j,k) fsdept_0(i,j,k)121 # define fsdepw_a(i,j,k) fsdepw_0(i,j,k)122 # define fsde3w_a(i,j,k) fsde3w_0(i,j,k)123 103 # define fse3t_a(i,j,k) fse3t_0(i,j,k) 124 104 # define fse3u_a(i,j,k) fse3u_0(i,j,k) 125 105 # define fse3v_a(i,j,k) fse3v_0(i,j,k) 126 # define fse3f_a(i,j,k) fse3f_0(i,j,k)127 # define fse3w_a(i,j,k) fse3w_0(i,j,k)128 # define fse3uw_a(i,j,k) fse3uw_0(i,j,k)129 # define fse3vw_a(i,j,k) fse3vw_0(i,j,k)130 106 #endif 131 107 !!---------------------------------------------------------------------- -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/divcur.F90
r1953 r2148 14 14 USE in_out_manager ! I/O manager 15 15 USE obc_oce ! ocean lateral open boundary condition 16 USE bdy_oce 16 USE bdy_oce ! Unstructured open boundaries variables 17 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 18 18 … … 87 87 INTEGER :: ii, ij, jl ! temporary integer 88 88 INTEGER :: ijt, iju ! temporary integer 89 REAL(wp) :: zraur, zdep 89 90 REAL(wp), DIMENSION( jpi ,1:jpj+2) :: zwu ! workspace 90 91 REAL(wp), DIMENSION(-1:jpi+2, jpj ) :: zwv ! workspace … … 301 302 !! * Local declarations 302 303 INTEGER :: ji, jj, jk ! dummy loop indices 304 REAL(wp) :: zraur, zdep 303 305 !!---------------------------------------------------------------------- 304 306 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynnxt.F90
r1970 r2148 22 22 USE oce ! ocean dynamics and tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE sbc_oce ! Surface boundary condition: ocean fields 25 USE phycst ! physical constants 24 26 USE dynspg_oce ! type of surface pressure gradient 25 27 USE dynadv ! dynamics: vector invariant versus flux form … … 87 89 !! un,vn now horizontal velocity of next time-step 88 90 !!---------------------------------------------------------------------- 91 USE oce, ONLY : ze3u_f => ta ! use ta as 3D workspace 92 USE oce, ONLY : ze3v_f => sa ! use sa as 3D workspace 89 93 INTEGER, INTENT( in ) :: kt ! ocean time-step index 90 94 !! … … 95 99 REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar 96 100 REAL(wp) :: zve3a , zve3n , zve3b ! - - 97 REAL(wp) :: ze3u_b, ze3u_n, ze3u_a ! - -98 REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - -99 101 REAL(wp) :: zuf , zvf ! - - 102 REAL(wp) :: zec ! - - 103 REAL(wp) :: zv_t_ij , zv_t_ip1j ! - - 104 REAL(wp) :: zv_t_ijp1 ! - - 105 REAL(wp), DIMENSION(jpi,jpj) :: zs_t, zs_u_1, zs_v_1 ! temporary 2D workspace 100 106 !!---------------------------------------------------------------------- 101 107 … … 146 152 # if defined key_obc 147 153 ! !* OBC open boundaries 148 IF( lk_obc )CALL obc_dyn( kt )154 CALL obc_dyn( kt ) 149 155 ! 150 156 IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN … … 212 218 END DO 213 219 ELSE !* Leap-Frog : Asselin filter and swap 214 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 220 ! ! =============! 221 IF( .NOT. lk_vvl ) THEN ! Fixed volume ! 222 ! ! =============! 215 223 DO jk = 1, jpkm1 216 224 DO jj = 1, jpj 217 225 DO ji = 1, jpi 218 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)219 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)226 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 227 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 220 228 ! 221 229 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 226 234 END DO 227 235 END DO 228 ELSE ! applied on thickness weighted velocity 236 ! ! ================! 237 ELSE ! Variable volume ! 238 ! ! ================! 239 ! Before scale factor at t-points 240 ! ------------------------------- 229 241 DO jk = 1, jpkm1 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ze3u_a = fse3u_a(ji,jj,jk) 233 ze3v_a = fse3v_a(ji,jj,jk) 234 ze3u_n = fse3u_n(ji,jj,jk) 235 ze3v_n = fse3v_n(ji,jj,jk) 236 ze3u_b = fse3u_b(ji,jj,jk) 237 ze3v_b = fse3v_b(ji,jj,jk) 238 ! 239 zue3a = ua(ji,jj,jk) * ze3u_a 240 zve3a = va(ji,jj,jk) * ze3v_a 241 zue3n = un(ji,jj,jk) * ze3u_n 242 zve3n = vn(ji,jj,jk) * ze3v_n 243 zue3b = ub(ji,jj,jk) * ze3u_b 244 zve3b = vb(ji,jj,jk) * ze3v_b 245 ! 246 zuf = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 247 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 248 zvf = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 249 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 250 ! 251 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 252 vb(ji,jj,jk) = zvf 253 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 254 vn(ji,jj,jk) = va(ji,jj,jk) 255 END DO 256 END DO 257 END DO 242 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) & 243 & + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) & 244 & - 2.e0 * fse3t_n(:,:,jk) ) 245 ENDDO 246 ! Add volume filter correction only at the first level of t-point scale factors 247 zec = atfp * rdt / rau0 248 fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 249 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 250 zs_t (:,:) = e1t(:,:) * e2t(:,:) 251 zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 252 zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 253 ! 254 IF( ln_dynadv_vec ) THEN 255 ! Before scale factor at (u/v)-points 256 ! ----------------------------------- 257 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 258 DO jk = 1, jpkm1 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 261 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 262 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 263 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 264 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 265 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 266 END DO 267 END DO 268 END DO 269 ! lateral boundary conditions 270 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 271 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 272 ! Add initial scale factor to scale factor anomaly 273 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 274 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 275 ! Leap-Frog - Asselin filter and swap: applied on velocity 276 ! ----------------------------------- 277 DO jk = 1, jpkm1 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 281 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 282 ! 283 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 284 vb(ji,jj,jk) = zvf 285 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 286 vn(ji,jj,jk) = va(ji,jj,jk) 287 END DO 288 END DO 289 END DO 290 ! 291 ELSE 292 ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 293 !----------------------------------------------- 294 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 295 DO jk = 1, jpkm1 296 DO jj = 1, jpjm1 297 DO ji = 1, jpim1 298 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 299 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 300 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 301 ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 302 ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 303 END DO 304 END DO 305 END DO 306 ! lateral boundary conditions 307 CALL lbc_lnk( ze3u_f, 'U', 1. ) 308 CALL lbc_lnk( ze3v_f, 'V', 1. ) 309 ! Add initial scale factor to scale factor anomaly 310 ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 311 ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 312 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 313 ! ----------------------------------- =========================== 314 DO jk = 1, jpkm1 315 DO jj = 1, jpj 316 DO ji = 1, jpim1 317 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 318 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 319 zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk) 320 zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk) 321 zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk) 322 zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 323 ! 324 zuf = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 325 zvf = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 326 ! 327 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 328 vb(ji,jj,jk) = zvf 329 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 330 vn(ji,jj,jk) = va(ji,jj,jk) 331 END DO 332 END DO 333 END DO 334 fse3u_b(:,:,:) = ze3u_f(:,:,:) ! e3u_b <-- filtered scale factor 335 fse3v_b(:,:,:) = ze3v_f(:,:,:) 336 CALL lbc_lnk( ub, 'U', -1. ) ! lateral boundary conditions 337 CALL lbc_lnk( vb, 'V', -1. ) 338 ENDIF 339 ! 258 340 ENDIF 341 ! 259 342 ENDIF 260 343 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r1537 r2148 4 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 5 5 !!============================================================================== 6 !! History : ! 90-10 (B. Blanke) Original code 7 !! ! 97-05 (G. Madec) vertical component of isopycnal 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal 8 !! NEMO 1.0 ! 1002-08 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! sion using an explicit time-stepping scheme. 14 15 !!---------------------------------------------------------------------- 15 !! * Modules used16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain … … 24 24 PRIVATE 25 25 26 !! * Routine accessibility 27 PUBLIC dyn_zdf_exp ! called by step.F90 26 PUBLIC dyn_zdf_exp ! called by step.F90 28 27 29 28 !! * Substitutions … … 31 30 # include "vectopt_loop_substitute.h90" 32 31 !!---------------------------------------------------------------------- 33 !! OPA 9.0 , LOCEAN-IPSL (2005)32 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 34 33 !! $Id$ 35 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 47 46 !! technique. The vertical diffusion of momentum is given by: 48 47 !! diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) ) 49 !! Surface boundary conditions: wind stress input 48 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 50 49 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F90) 51 50 !! Add this trend to the general trend ua : … … 54 53 !! ** Action : - Update (ua,va) with the vertical diffusive trend 55 54 !!--------------------------------------------------------------------- 56 !! * Arguments 57 INTEGER , INTENT( in ) :: kt ! ocean time-step index 58 REAL(wp), INTENT( in ) :: p2dt ! time-step 59 60 !! * Local declarations 61 INTEGER :: ji, jj, jk, jl ! dummy loop indices 62 REAL(wp) :: zrau0r, zlavmr, zua, zva ! temporary scalars 63 REAL(wp), DIMENSION(jpi,jpk) :: zwx, zwy, zwz, zww ! temporary workspace arrays 55 INTEGER , INTENT(in) :: kt ! ocean time-step index 56 REAL(wp), INTENT(in) :: p2dt ! time-step 57 !! 58 INTEGER :: ji, jj, jk, jl ! dummy loop indices 59 REAL(wp) :: zrau0r, zlavmr, zua, zva ! temporary scalars 60 REAL(wp), DIMENSION(jpi,jpk) :: zwx, zwy, zwz, zww ! 2D workspace 64 61 !!---------------------------------------------------------------------- 65 62 66 63 IF( kt == nit000 ) THEN 67 64 IF(lwp) WRITE(numout,*) 68 IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'65 IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator' 69 66 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 70 67 ENDIF 71 68 72 ! Local constant initialization 73 ! ----------------------------- 74 zrau0r = 1. / rau0 ! inverse of the reference density 75 zlavmr = 1. / float( nn_zdfexp ) ! inverse of the number of sub time step 69 zrau0r = 1. / rau0 ! Local constant initialization 70 zlavmr = 1. / REAL( nn_zdfexp ) 76 71 77 ! ! =============== 78 DO jj = 2, jpjm1 ! Vertical slab 79 ! ! =============== 80 81 ! Surface boundary condition 82 DO ji = 2, jpim1 83 zwy(ji,1) = utau(ji,jj) * zrau0r 84 zww(ji,1) = vtau(ji,jj) * zrau0r 72 ! ! =============== 73 DO jj = 2, jpjm1 ! Vertical slab 74 ! ! =============== 75 DO ji = 2, jpim1 ! Surface boundary condition 76 zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r 77 zww(ji,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r 85 78 END DO 86 87 ! Initialization of x, z and contingently trends array 88 DO jk = 1, jpk 79 DO jk = 1, jpk ! Initialization of x, z and contingently trends array 89 80 DO ji = 2, jpim1 90 81 zwx(ji,jk) = ub(ji,jj,jk) … … 92 83 END DO 93 84 END DO 94 95 ! Time splitting loop 96 DO jl = 1, nn_zdfexp 97 98 ! First vertical derivative 99 DO jk = 2, jpk 85 ! 86 DO jl = 1, nn_zdfexp ! Time splitting loop 87 ! 88 DO jk = 2, jpk ! First vertical derivative 100 89 DO ji = 2, jpim1 101 90 zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) … … 103 92 END DO 104 93 END DO 105 106 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 107 DO jk = 1, jpkm1 94 DO jk = 1, jpkm1 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 108 95 DO ji = 2, jpim1 109 zua = zlavmr *( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)110 zva = zlavmr *( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)96 zua = zlavmr * ( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk) 97 zva = zlavmr * ( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk) 111 98 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 112 99 va(ji,jj,jk) = va(ji,jj,jk) + zva 113 114 zwx(ji,jk) = zwx(ji,jk) + p2dt *zua*umask(ji,jj,jk)115 zwz(ji,jk) = zwz(ji,jk) + p2dt *zva*vmask(ji,jj,jk)100 ! 101 zwx(ji,jk) = zwx(ji,jk) + p2dt * zua * umask(ji,jj,jk) 102 zwz(ji,jk) = zwz(ji,jk) + p2dt * zva * vmask(ji,jj,jk) 116 103 END DO 117 104 END DO 118 105 ! 119 106 END DO 120 121 ! ! =============== 122 END DO ! End of slab 123 ! ! =============== 124 107 ! ! =============== 108 END DO ! End of slab 109 ! ! =============== 125 110 END SUBROUTINE dyn_zdf_exp 126 111 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r1662 r2148 4 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 5 5 !!============================================================================== 6 !! History : ! 90-10 (B. Blanke) Original code 7 !! ! 97-05 (G. Madec) vertical component of isopycnal 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal 8 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! sion using a implicit time-stepping. 14 15 !!---------------------------------------------------------------------- 15 !! OPA 9.0 , LOCEAN-IPSL (2005)16 !! $Id$17 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)18 !!----------------------------------------------------------------------19 !! * Modules used20 16 USE oce ! ocean dynamics and tracers 21 17 USE dom_oce ! ocean space and time domain … … 28 24 PRIVATE 29 25 30 !! * Routine accessibility 31 PUBLIC dyn_zdf_imp ! called by step.F90 26 PUBLIC dyn_zdf_imp ! called by step.F90 32 27 33 28 !! * Substitutions … … 35 30 # include "vectopt_loop_substitute.h90" 36 31 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (2005)32 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 38 33 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 35 !!---------------------------------------------------------------------- 41 36 42 37 CONTAINS 43 44 38 45 39 SUBROUTINE dyn_zdf_imp( kt, p2dt ) … … 54 48 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 55 49 !! backward time stepping 56 !! Surface boundary conditions: wind stress input 50 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 57 51 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 58 52 !! Add this trend to the general trend ua : 59 53 !! ua = ua + dz( avmu dz(u) ) 60 54 !! 61 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 62 !! mixing trend. 55 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 63 56 !!--------------------------------------------------------------------- 64 !! * Modules used 65 USE oce, ONLY : zwd => ta, & ! use ta as workspace 66 zws => sa ! use sa as workspace 67 68 !! * Arguments 69 INTEGER , INTENT( in ) :: kt ! ocean time-step index 70 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 71 72 !! * Local declarations 73 INTEGER :: ji, jj, jk ! dummy loop indices 74 REAL(wp) :: zrau0r, z2dtf, zcoef, zzws, zrhs ! temporary scalars 75 REAL(wp) :: zzwi ! temporary scalars 76 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! temporary workspace arrays 57 USE oce, ONLY : zwd => ta ! use ta as workspace 58 USE oce, ONLY : zws => sa ! use sa as workspace 59 !! 60 INTEGER , INTENT( in ) :: kt ! ocean time-step index 61 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 62 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: zrau0r, zcoef ! temporary scalars 65 REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace 77 67 !!---------------------------------------------------------------------- 78 68 … … 95 85 ! friction velocity in dyn_bfr using values from the previous timestep. There 96 86 ! is no need to include these in the implicit calculation. 97 DO jk = 1, jpkm1 87 ! 88 DO jk = 1, jpkm1 ! Matrix 98 89 DO jj = 2, jpjm1 99 90 DO ji = fs_2, fs_jpim1 ! vector opt. 100 91 zcoef = - p2dt / fse3u(ji,jj,jk) 101 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )102 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)103 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)104 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)92 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 93 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 94 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 95 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 105 96 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 106 97 END DO 107 98 END DO 108 99 END DO 109 110 ! Surface boudary conditions 111 DO jj = 2, jpjm1 100 DO jj = 2, jpjm1 ! Surface boudary conditions 112 101 DO ji = fs_2, fs_jpim1 ! vector opt. 113 102 zwi(ji,jj,1) = 0. … … 130 119 ! The solution (the after velocity) is in ua 131 120 !----------------------------------------------------------------------- 132 133 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 134 DO jk = 2, jpkm1 121 ! 122 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 135 123 DO jj = 2, jpjm1 136 124 DO ji = fs_2, fs_jpim1 ! vector opt. … … 139 127 END DO 140 128 END DO 141 142 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 143 DO jj = 2, jpjm1 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 146 !!! ua(ji,jj,1) = ub(ji,jj,1) & 147 !!! + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) ) 148 z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 ) 149 ua(ji,jj,1) = ub(ji,jj,1) & 150 + p2dt * ua(ji,jj,1) + z2dtf * utau(ji,jj) 129 ! 130 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 131 DO ji = fs_2, fs_jpim1 ! vector opt. 132 ua(ji,jj,1) = ub(ji,jj,1) & 133 & + p2dt * ( ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) ) 151 134 END DO 152 135 END DO … … 159 142 END DO 160 143 END DO 161 162 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 163 DO jj = 2, jpjm1 144 ! 145 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 164 146 DO ji = fs_2, fs_jpim1 ! vector opt. 165 147 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 173 155 END DO 174 156 END DO 175 176 ! Normalization to obtain the general momentum trend ua 177 DO jk = 1, jpkm1 157 ! 158 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend ua == 178 159 DO jj = 2, jpjm1 179 160 DO ji = fs_2, fs_jpim1 ! vector opt. … … 192 173 ! friction velocity in dyn_bfr using values from the previous timestep. There 193 174 ! is no need to include these in the implicit calculation. 194 DO jk = 1, jpkm1 175 ! 176 DO jk = 1, jpkm1 ! Matrix 195 177 DO jj = 2, jpjm1 196 178 DO ji = fs_2, fs_jpim1 ! vector opt. 197 179 zcoef = -p2dt / fse3v(ji,jj,jk) 198 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )180 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 199 181 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 200 zzws = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)182 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 201 183 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 202 184 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws … … 204 186 END DO 205 187 END DO 206 207 ! Surface boudary conditions 208 DO jj = 2, jpjm1 188 DO jj = 2, jpjm1 ! Surface boudary conditions 209 189 DO ji = fs_2, fs_jpim1 ! vector opt. 210 190 zwi(ji,jj,1) = 0.e0 … … 223 203 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 224 204 ! 225 ! m is decomposed in the product of an upper and lower triangular 226 ! matrix 205 ! m is decomposed in the product of an upper and lower triangular matrix 227 206 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 228 207 ! The solution (after velocity) is in 2d array va 229 208 !----------------------------------------------------------------------- 230 231 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 232 DO jk = 2, jpkm1 209 ! 210 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 233 211 DO jj = 2, jpjm1 234 212 DO ji = fs_2, fs_jpim1 ! vector opt. … … 237 215 END DO 238 216 END DO 239 240 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 244 !!! va(ji,jj,1) = vb(ji,jj,1) & 245 !!! + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) ) 246 z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 ) 247 va(ji,jj,1) = vb(ji,jj,1) & 248 + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj) 217 ! 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 va(ji,jj,1) = vb(ji,jj,1) & 221 & + p2dt * ( va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) ) 249 222 END DO 250 223 END DO … … 257 230 END DO 258 231 END DO 259 260 ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 261 DO jj = 2, jpjm1 232 ! 233 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 262 234 DO ji = fs_2, fs_jpim1 ! vector opt. 263 235 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 271 243 END DO 272 244 END DO 273 274 ! Normalization to obtain the general momentum trend va 275 DO jk = 1, jpkm1 245 ! 246 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend va == 276 247 DO jj = 2, jpjm1 277 248 DO ji = fs_2, fs_jpim1 ! vector opt. … … 280 251 END DO 281 252 END DO 282 253 ! 283 254 END SUBROUTINE dyn_zdf_imp 284 255 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90
r2113 r2148 5 5 !!============================================================================== 6 6 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 7 8 !!---------------------------------------------------------------------- 8 9 … … 27 28 USE diaar5, ONLY : lk_diaar5 28 29 USE iom 29 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 31 31 32 IMPLICIT NONE … … 38 39 # include "domzgr_substitute.h90" 39 40 # include "vectopt_loop_substitute.h90" 40 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 43 43 !! $Id$ 44 44 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 54 54 !! and update the now vertical coordinate (lk_vvl=T). 55 55 !! 56 !! ** Method : - 57 !! 58 !! - Using the incompressibility hypothesis, the vertical 56 !! ** Method : - Using the incompressibility hypothesis, the vertical 59 57 !! velocity is computed by integrating the horizontal divergence 60 58 !! from the bottom to the surface minus the scale factor evolution. … … 63 61 !! ** action : ssha : after sea surface height 64 62 !! wn : now vertical velocity 65 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 66 !! at u-, v-, f-point s 67 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 63 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T) 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !! 66 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 68 67 !!---------------------------------------------------------------------- 69 68 USE oce, ONLY : z3d => ta ! use ta as 3D workspace … … 73 72 INTEGER :: ji, jj, jk ! dummy loop indices 74 73 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 75 REAL(wp) :: z2dt, z raur! temporary scalars74 REAL(wp) :: z2dt, z1_2dt, z1_rau0 ! temporary scalars 76 75 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 77 76 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 79 78 80 79 IF( kt == nit000 ) THEN 80 ! 81 81 IF(lwp) WRITE(numout,*) 82 82 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' … … 95 95 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 96 96 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 97 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &98 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) )99 97 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 100 98 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 101 99 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 102 100 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 103 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &104 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )105 101 END DO 106 102 END DO 107 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 108 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 109 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 105 DO jj = 1, jpjm1 106 DO ji = 1, jpim1 ! NO Vector Opt. 107 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 108 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 109 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 110 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 111 END DO 112 END DO 113 CALL lbc_lnk( sshf_n, 'F', 1. ) 110 114 ENDIF 111 115 ! 112 116 ENDIF 113 117 114 ! !------------------------------ !115 IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)116 ! !------------------------------!118 ! !------------------------------------------! 119 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case) 120 ! !------------------------------------------! 117 121 DO jk = 1, jpkm1 118 fsdept(:,:,jk) = fsdept_n(:,:,jk) 122 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays 119 123 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 120 124 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 121 125 ! 122 fse3t (:,:,jk) = fse3t_n (:,:,jk) 126 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays 123 127 fse3u (:,:,jk) = fse3u_n (:,:,jk) 124 128 fse3v (:,:,jk) = fse3v_n (:,:,jk) … … 128 132 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 129 133 END DO 130 ! ! now ocean depth (at u- and v-points)131 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 134 ! 135 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 132 136 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 133 ! 137 ! ! now masked inverse of the ocean depth (at u- and v-points) 134 138 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 135 139 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 136 140 ! 137 ! 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 rnf_dep(ji,jj) = 0. 141 DO jk = 1, rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth 142 rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box 143 ENDDO 144 ENDDO 145 ENDDO 146 ! 147 ENDIF 148 149 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 150 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 151 152 ! set time step size (Euler/Leapfrog) 153 z2dt = 2. * rdt 141 ENDIF 142 ! 143 144 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 145 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 146 147 z2dt = 2. * rdt ! set time step size (Euler/Leapfrog) 154 148 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 155 156 zraur = 1. / rau0157 149 158 150 ! !------------------------------! … … 163 155 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 164 156 END DO 165 166 157 ! ! Sea surface elevation time stepping 167 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 158 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 159 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b -2*rnf ) = emp 160 z1_rau0 = 0.5 / rau0 161 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) - 2 * rnf(:,:) ) + zhdiv(:,:) ) ) & 162 & * tmask(:,:,1) 168 163 169 164 #if defined key_obc 170 IF 165 IF( Agrif_Root() ) THEN 171 166 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 172 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)167 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 173 168 ENDIF 174 169 #endif 175 176 170 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 177 171 IF( lk_vvl ) THEN ! (required only in key_vvl case) … … 184 178 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 185 179 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 186 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &187 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) &188 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) )189 180 END DO 190 181 END DO 191 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 182 ! Boundaries conditions 183 CALL lbc_lnk( sshu_a, 'U', 1. ) 192 184 CALL lbc_lnk( sshv_a, 'V', 1. ) 193 CALL lbc_lnk( sshf_a, 'F', 1. ) 194 ENDIF 195 185 ENDIF 196 186 ! !------------------------------! 197 187 ! ! Now Vertical Velocity ! 198 188 ! !------------------------------! 199 ! ! integrate from the bottom the hor. divergence 200 DO jk = jpkm1, 1, -1 201 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 202 & - ( fse3t_a(:,:,jk) & 203 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 189 z1_2dt = 1.e0 / z2dt 190 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 192 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 193 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 194 & * tmask(:,:,jk) * z1_2dt 204 195 END DO 205 ! 196 197 ! !------------------------------! 198 ! ! outputs ! 199 ! !------------------------------! 206 200 CALL iom_put( "woce", wn ) ! vertical velocity 207 201 CALL iom_put( "ssh" , sshn ) ! sea surface height 208 202 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 209 IF( lk_diaar5 ) THEN 203 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 204 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 210 205 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 211 206 DO jk = 1, jpk 212 207 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 213 208 END DO 214 CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport 215 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport 216 ENDIF 209 CALL iom_put( "w_masstr" , z3d ) 210 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 211 ENDIF 212 ! 213 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 217 214 ! 218 215 END SUBROUTINE ssh_wzv … … 227 224 !! ssha already computed in ssh_wzv 228 225 !! 229 !! ** Method : - apply Asselin time fiter to now ssh and swap : 230 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 231 !! sshn = ssha 226 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 227 !! from the filter, see Leclair and Madec 2010) and swap : 228 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 229 !! - atfp * rdt * ( emp_b - emp ) / rau0 230 !! sshn = ssha 232 231 !! 233 232 !! ** action : - sshb, sshn : before & now sea surface height 234 233 !! ready for the next time step 235 !!---------------------------------------------------------------------- 236 INTEGER, INTENT( in ) :: kt ! ocean time-step index 237 !! 238 INTEGER :: ji, jj ! dummy loop indices 234 !! 235 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 236 !!---------------------------------------------------------------------- 237 INTEGER, INTENT(in) :: kt ! ocean time-step index 238 !! 239 INTEGER :: ji, jj ! dummy loop indices 240 REAL(wp) :: zec ! temporary scalar 239 241 !!---------------------------------------------------------------------- 240 242 … … 245 247 ENDIF 246 248 247 ! Time filter and swap of the ssh 248 ! ------------------------------- 249 250 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 251 ! ! ---------------------- ! 249 ! !--------------------------! 250 IF( lk_vvl ) THEN ! Variable volume levels ! 251 ! !--------------------------! 252 ! 253 ! ssh at t-, u-, v, f-points 254 !=========================== 252 255 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 253 256 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 254 257 sshu_n(:,:) = sshu_a(:,:) 255 258 sshv_n(:,:) = sshv_a(:,:) 256 sshf_n(:,:) = sshf_a(:,:) 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 ! NO Vector Opt. 261 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 262 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 263 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 264 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 265 END DO 266 END DO 267 ! Boundaries conditions 268 CALL lbc_lnk( sshf_n, 'F', 1. ) 257 269 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 270 zec = atfp * rdt / rau0 258 271 DO jj = 1, jpj 259 272 DO ji = 1, jpi ! before <-- now filtered 260 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 261 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 262 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 263 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 273 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 274 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 264 275 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 265 276 sshu_n(ji,jj) = sshu_a(ji,jj) 266 277 sshv_n(ji,jj) = sshv_a(ji,jj) 267 sshf_n(ji,jj) = sshf_a(ji,jj) 268 END DO 269 END DO 278 END DO 279 END DO 280 DO jj = 1, jpjm1 281 DO ji = 1, jpim1 ! NO Vector Opt. 282 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 283 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 284 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 285 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 286 END DO 287 END DO 288 ! Boundaries conditions 289 CALL lbc_lnk( sshf_n, 'F', 1. ) 290 DO jj = 1, jpjm1 291 DO ji = 1, jpim1 ! NO Vector Opt. 292 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 293 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 294 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 295 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 296 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 297 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 298 END DO 299 END DO 300 ! Boundaries conditions 301 CALL lbc_lnk( sshu_b, 'U', 1. ) 302 CALL lbc_lnk( sshv_b, 'V', 1. ) 270 303 ENDIF 271 ! 272 ELSE ! fixed levels : ssh at t-point only 273 ! ! ------------ ! 304 ! !--------------------------! 305 ELSE ! fixed levels ! 306 ! !--------------------------! 307 ! 308 ! ssh at t-point only 309 !==================== 274 310 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 275 311 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) … … 278 314 DO jj = 1, jpj 279 315 DO ji = 1, jpi ! before <-- now filtered 280 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 316 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 281 317 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 282 318 END DO … … 286 322 ENDIF 287 323 ! 288 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )324 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 289 325 ! 290 326 END SUBROUTINE ssh_nxt -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/IOM/restart.F90
r2104 r2148 7 7 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form 8 8 !! 2.0 ! 2006-07 (S. Masson) use IOM for restart 9 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D) 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 10 !! - - ! 2010-10 (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D) 10 11 !!---------------------------------------------------------------------- 11 12 … … 26 27 USE zdfmxl ! mixed layer depth 27 28 USE trdmld_oce ! ocean active mixed layer tracers trends variables 29 USE domvvl ! variable volume 28 30 USE traswp ! swap from 4D T-S to 3D T & S and vice versa 29 31 … … 39 41 40 42 !! * Substitutions 43 # include "domzgr_substitute.h90" 41 44 # include "vectopt_loop_substitute.h90" 42 45 !!---------------------------------------------------------------------- … … 123 126 CALL iom_rstput( kt, nitrst, numrow, 'hdivb' , hdivb ) 124 127 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb ) 128 IF( lk_vvl ) & 129 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 125 130 ! 126 131 CALL iom_rstput( kt, nitrst, numrow, 'un' , un ) ! now fields … … 154 159 !!---------------------------------------------------------------------- 155 160 REAL(wp) :: zrdt, zrdttra1 156 INTEGER :: j libalt = jprstlib161 INTEGER :: jk, jlibalt = jprstlib 157 162 LOGICAL :: llok 158 163 !!---------------------------------------------------------------------- … … 192 197 CALL iom_get( numror, jpdom_autoglo, 'hdivb', hdivb ) 193 198 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb ) 199 IF( lk_vvl ) & 200 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 194 201 ! 195 202 CALL iom_get( numror, jpdom_autoglo, 'un' , un ) ! now fields … … 219 226 hdivb(:,:,:) = hdivn(:,:,:) 220 227 sshb (:,:) = sshn (:,:) 228 IF( lk_vvl ) THEN 229 DO jk = 1, jpk 230 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 231 ENDDO 232 ENDIF 221 233 ENDIF 222 234 ! -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbc_oce.F90
r2000 r2148 6 6 !! History : 3.0 ! 2006-06 (G. Madec) Original code 7 7 !! - ! 2008-08 (G. Madec) namsbc moved from sbcmod 8 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 8 9 !!---------------------------------------------------------------------- 9 10 USE par_oce ! ocean parameters … … 37 38 !! Ocean Surface Boundary Condition fields 38 39 !!---------------------------------------------------------------------- 39 LOGICAL , PUBLIC :: lhftau = .FALSE. !: HF tau contribution: mean of stress module - module of the mean stress 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: utau !: sea surface i-stress (ocean referential) [N/m2] 41 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: vtau !: sea surface j-stress (ocean referential) [N/m2] 42 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: taum !: module of sea surface stress (at T-point) [N/m2] 43 !! wndm is used only in PISCES to compute gases exchanges at the surface of the free ocean or in the leads in sea-ice parts 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s] 45 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qsr !: sea heat flux: solar [W/m2] 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qns !: sea heat flux: non solar [W/m2] 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2] 48 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2] 49 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emp !: freshwater budget: volume flux [Kg/m2/s] 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emps !: freshwater budget: concentration/dillution [Kg/m2/s] 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rnf !: river runoff [Kg/m2/s] 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emp_tot !: total evaporation - (liquid + solid) precpitation over oce and ice 53 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: tprecip !: total precipitation [Kg/m2/s] 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sprecip !: solid precipitation [Kg/m2/s] 55 !!$ REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rrunoff !: runoff 56 !!$ REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: calving !: calving 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fr_i !: ice fraction (between 0 to 1) - 40 LOGICAL , PUBLIC :: lhftau = .FALSE. !: HF tau used in TKE: mean(stress module) - module(mean stress) 41 !! !! now ! before !! 42 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: utau , utau_b !: sea surface i-stress (ocean referential) [N/m2] 43 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: vtau , vtau_b !: sea surface j-stress (ocean referential) [N/m2] 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: taum !: module of sea surface stress (at T-point) [N/m2] 45 !! wndm is used only in PISCES to compute surface gases exchanges in ice-free ocean or leads 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s] 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qsr !: sea heat flux: solar [W/m2] 48 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qns , qns_b !: sea heat flux: non solar [W/m2] 49 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2] 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2] 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s] 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emps , emps_b !: freshwater budget: concentration/dillution [Kg/m2/s] 53 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s] 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rnf !: river runoff [Kg/m2/s] 55 ! - ML - begin 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_hc_n !: sbc heat content trend now [K.m/s] 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_hc_b !: " " " " before " 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_sc_n !: sbc salt content trend now [psu.m/s] 59 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_sc_b !: " " " " before " 60 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: qsr_hc_n !: heat content trend due to qsr flux now [K.m/s] 61 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: qsr_hc_b !: " " " " " " " before " 62 ! - ML - end 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: tprecip !: total precipitation [Kg/m2/s] 64 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sprecip !: solid precipitation [Kg/m2/s] 65 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1) 58 66 #if defined key_cpl_carbon_cycle 59 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: atm_co2 !: atmospheric pCO2 [ppm]67 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: atm_co2 !: atmospheric pCO2 [ppm] 60 68 #endif 69 !!$ REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rrunoff !: runoff 70 !!$ REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: calving !: calving 61 71 62 72 !!---------------------------------------------------------------------- … … 71 81 72 82 !!---------------------------------------------------------------------- 73 !! OPA 9.0 , LOCEAN-IPSL (2006)74 !! $ Id:$83 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 84 !! $Id$ 75 85 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 76 86 !!====================================================================== 77 78 87 END MODULE sbc_oce -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcmod.F90
r2000 r2148 4 4 !! Surface module : provide to the ocean its surface boundary condition 5 5 !!====================================================================== 6 !! History : 3.0 ! 07-2006 (G. Madec) Original code 7 !! - ! 08-2008 (S. Masson, E. .... ) coupled interface 6 !! History : 3.0 ! 2006-07 (G. Madec) Original code 7 !! 3.1 ! 2008-08 (S. Masson, E. Maisonnave, G. Madec) coupled interface 8 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 8 9 !!---------------------------------------------------------------------- 9 10 … … 49 50 # include "domzgr_substitute.h90" 50 51 !!---------------------------------------------------------------------- 51 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)52 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 52 53 !! $Id$ 53 54 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 86 87 !!gm IF( lk_sbc_cpl ) THEN ; ln_cpl = .TRUE. ; ELSE ; ln_cpl = .FALSE. ; ENDIF 87 88 88 IF 89 IF( Agrif_Root() ) THEN 89 90 IF( lk_lim2 ) nn_ice = 2 90 91 IF( lk_lim3 ) nn_ice = 3 … … 179 180 !! CAUTION : never mask the surface stress field (tke sbc) 180 181 !! 181 !! ** Action : - set the ocean surface boundary condition, i.e. 182 !! utau, vtau, qns, qsr, emp, emps, qrp, erp 182 !! ** Action : - set the ocean surface boundary condition at before and now 183 !! time step, i.e. 184 !! utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b 185 !! utau , vtau , qns , qsr , emp , emps , qrp , erp 183 186 !! - updte the ice fraction : fr_i 184 187 !!---------------------------------------------------------------------- … … 186 189 !!--------------------------------------------------------------------- 187 190 188 CALL iom_setkt( kt + nn_fsbc - 1 ) ! in sbc, iom_put is called every nn_fsbc time step 189 ! 190 ! ocean to sbc mean sea surface variables (ss._m) 191 ! --------------------------------------- 192 CALL sbc_ssm( kt ) ! sea surface mean currents (at U- and V-points), 193 ! ! temperature and salinity (at T-point) over nf_sbc time-step 194 ! ! (i.e. sst_m, sss_m, ssu_m, ssv_m) 195 196 ! sbc formulation 197 ! --------------- 198 199 SELECT CASE( nsbc ) ! Compute ocean surface boundary condition 200 ! ! (i.e. utau,vtau, qns, qsr, emp, emps) 191 ! ! ---------------------------------------- ! 192 IF( kt /= nit000 ) THEN ! Swap of forcing fields ! 193 ! ! ---------------------------------------- ! 194 utau_b(:,:) = utau(:,:) ! Swap the ocean forcing fields 195 vtau_b(:,:) = vtau(:,:) ! (except at nit000 where before fields 196 qns_b (:,:) = qns (:,:) ! are set at the end of the routine) 197 ! The 3D heat content due to qsr forcing is treated in traqsr 198 ! qsr_b (:,:) = qsr (:,:) 199 emp_b (:,:) = emp (:,:) 200 emps_b(:,:) = emps(:,:) 201 ENDIF 202 ! ! ---------------------------------------- ! 203 ! ! forcing field computation ! 204 ! ! ---------------------------------------- ! 205 206 CALL iom_setkt( kt + nn_fsbc - 1 ) ! in sbc, iom_put is called every nn_fsbc time step 207 ! 208 CALL sbc_ssm( kt ) ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 209 ! ! averaged over nf_sbc time-step 210 211 !== sbc formulation ==! 212 213 SELECT CASE( nsbc ) ! Compute ocean surface boundary condition 214 ! ! (i.e. utau,vtau, qns, qsr, emp, emps) 201 215 CASE( 0 ) ; CALL sbc_gyre ( kt ) ! analytical formulation : GYRE configuration 202 216 CASE( 1 ) ; CALL sbc_ana ( kt ) ! analytical formulation : uniform sbc … … 214 228 END SELECT 215 229 216 ! Misc. Options 217 ! ------------- 230 ! !== Misc. Options ==! 218 231 219 232 !!gm IF( ln_dm2dc ) CALL sbc_dcy( kt ) ! Daily mean qsr distributed over the Diurnal Cycle … … 236 249 ! ! (update freshwater fluxes) 237 250 ! 251 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 252 ! ! ---------------------------------------- ! 253 IF( ln_rstart .AND. & !* Restart: read in restart file 254 & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 255 IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields red in the restart file' 256 CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b ) ! before i-stress (U-point) 257 CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b ) ! before j-stress (V-point) 258 CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b ) ! before non solar heat flux (T-point) 259 ! The 3D heat content due to qsr forcing is treated in traqsr 260 ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b ) ! before solar heat flux (T-point) 261 CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b ) ! before freshwater flux (T-point) 262 CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b ) ! before C/D freshwater flux (T-point) 263 ELSE !* no restart: set from nit000 values 264 IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields set to nit000' 265 utau_b(:,:) = utau(:,:) 266 vtau_b(:,:) = vtau(:,:) 267 qns_b (:,:) = qns (:,:) 268 ! qsr_b (:,:) = qsr (:,:) 269 emp_b (:,:) = emp (:,:) 270 emps_b(:,:) = emps(:,:) 271 ENDIF 272 ENDIF 273 ! ! ---------------------------------------- ! 274 IF( lrst_oce ) THEN ! Write in the ocean restart file ! 275 ! ! ---------------------------------------- ! 276 IF(lwp) WRITE(numout,*) 277 IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ', & 278 & 'at it= ', kt,' date= ', ndastp 279 IF(lwp) WRITE(numout,*) '~~~~' 280 CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau ) 281 CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau ) 282 CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns ) 283 ! The 3D heat content due to qsr forcing is treated in traqsr 284 ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr ) 285 CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp ) 286 CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps ) 287 ENDIF 288 289 ! ! ---------------------------------------- ! 290 ! ! Outputs and control print ! 291 ! ! ---------------------------------------- ! 238 292 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 239 CALL iom_put( "emp-rnf" , (emp-rnf) ) ! upward water flux240 CALL iom_put( "emps-rnf" , (emps-rnf) ) ! c/d water flux241 CALL iom_put( "qns+qsr" , qns + qsr )! total heat flux (caution if ln_dm2dc=true, to be242 CALL iom_put( "qns" , qns )! solar heat flux moved after the call to iom_setkt)243 CALL iom_put( "qsr" , qsr )! solar heat flux moved after the call to iom_setkt)244 IF( nn_ice > 0 ) CALL iom_put( "ice_cover", fr_i )! ice fraction293 CALL iom_put( "emp-rnf" , emp - rnf ) ! upward water flux 294 CALL iom_put( "emps-rnf", emps - rnf ) ! c/d water flux 295 CALL iom_put( "qns+qsr" , qns + qsr ) ! total heat flux (caution if ln_dm2dc=true, to be 296 CALL iom_put( "qns" , qns ) ! solar heat flux moved after the call to iom_setkt) 297 CALL iom_put( "qsr" , qsr ) ! solar heat flux moved after the call to iom_setkt) 298 IF( nn_ice > 0 ) CALL iom_put( "ice_cover", fr_i ) ! ice fraction 245 299 ENDIF 246 300 ! … … 253 307 ! 254 308 IF(ln_ctl) THEN ! print mean trends (used for debugging) 255 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i- : ', mask1=tmask, ovlap=1 )256 CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf - : ', mask1=tmask, ovlap=1 ) 257 CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 ) 258 CALL prt_ctl(tab2d_1=qns , clinfo1=' qns- : ', mask1=tmask, ovlap=1 )259 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr- : ', mask1=tmask, ovlap=1 )260 CALL prt_ctl(tab3d_1=tmask , clinfo1=' tmask: ', mask1=tmask, ovlap=1, kdim=jpk )261 CALL prt_ctl(tab3d_1=tn , clinfo1=' sst- : ', mask1=tmask, ovlap=1, kdim=1 )262 CALL prt_ctl(tab3d_1=sn , clinfo1=' sss- : ', mask1=tmask, ovlap=1, kdim=1 )263 CALL prt_ctl(tab2d_1=utau , clinfo1=' utau- : ', mask1=umask, &264 & tab2d_2=vtau , clinfo2=' vtau- : ', mask2=vmask, ovlap=1 )309 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 ) 310 CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf - : ', mask1=tmask, ovlap=1 ) 311 CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 ) 312 CALL prt_ctl(tab2d_1=qns , clinfo1=' qns - : ', mask1=tmask, ovlap=1 ) 313 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 ) 314 CALL prt_ctl(tab3d_1=tmask , clinfo1=' tmask - : ', mask1=tmask, ovlap=1, kdim=jpk ) 315 CALL prt_ctl(tab3d_1=tn , clinfo1=' sst - : ', mask1=tmask, ovlap=1, kdim=1 ) 316 CALL prt_ctl(tab3d_1=sn , clinfo1=' sss - : ', mask1=tmask, ovlap=1, kdim=1 ) 317 CALL prt_ctl(tab2d_1=utau , clinfo1=' utau - : ', mask1=umask, & 318 & tab2d_2=vtau , clinfo2=' vtau - : ', mask2=vmask, ovlap=1 ) 265 319 ENDIF 266 320 ! -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcrnf.F90
r2113 r2148 33 33 TYPE(FLD_N) , PUBLIC :: sn_cnf !: information about the runoff mouth file to be read 34 34 TYPE(FLD_N) :: sn_sal_rnf !: information about the salinities of runoff file to be read 35 TYPE(FLD_N) :: sn_t em_rnf !: information about the temperatures of runoff file to be read35 TYPE(FLD_N) :: sn_tmp_rnf !: information about the temperatures of runoff file to be read 36 36 TYPE(FLD_N) :: sn_dep_rnf !: information about the depth which river inflow affects 37 37 LOGICAL , PUBLIC :: ln_rnf_mouth = .false. !: specific treatment in mouths vicinity … … 53 53 INTEGER, PUBLIC, DIMENSION(jpi,jpj) :: rnf_mod_dep !: depth of runoff in model levels 54 54 REAL, PUBLIC, DIMENSION(jpi,jpj) :: rnf_sal !: salinity of river runoff 55 REAL, PUBLIC, DIMENSION(jpi,jpj) :: rnf_t em!: temperature of river runoff55 REAL, PUBLIC, DIMENSION(jpi,jpj) :: rnf_tmp !: temperature of river runoff 56 56 57 57 INTEGER :: ji, jj ,jk ! dummy loop indices … … 86 86 !!---------------------------------------------------------------------- 87 87 ! 88 IF( kt == nit000 ) THEN 89 CALL sbc_rnf_init ! Read namelist and allocate structures 90 ENDIF 88 IF( kt == nit000 ) CALL sbc_rnf_init ! Read namelist and allocate structures 91 89 92 90 ! !-------------------! … … 117 115 IF ( ln_rnf_att ) THEN 118 116 rnf_sal(:,:) = ( sf_sal_rnf(1)%fnow(:,:,1) ) 119 rnf_t em(:,:) = ( sf_tem_rnf(1)%fnow(:,:,1) )117 rnf_tmp(:,:) = ( sf_tem_rnf(1)%fnow(:,:,1) ) 120 118 ELSE 121 119 rnf_sal(:,:) = 0 122 rnf_t em(:,:) = -999120 rnf_tmp(:,:) = -999 123 121 ENDIF 124 122 CALL iom_put( "runoffs", rnf ) ! runoffs … … 143 141 CHARACTER(len=32) :: rn_dep_file ! runoff file name 144 142 !! 145 NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, sn_rnf, sn_cnf, sn_sal_rnf, sn_t em_rnf, sn_dep_rnf, &143 NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, sn_rnf, sn_cnf, sn_sal_rnf, sn_tmp_rnf, sn_dep_rnf, & 146 144 & ln_rnf_mouth, ln_rnf_att, rn_hrnf, rn_avt_rnf, rn_rfact 147 145 !!---------------------------------------------------------------------- … … 157 155 158 156 sn_sal_rnf = FLD_N( 'runoffs', 24. , 'rosaline' , .TRUE. , .true. , 'yearly' , '' , '' ) 159 sn_t em_rnf = FLD_N( 'runoffs', 24. , 'rotemper' , .TRUE. , .true. , 'yearly' , '' , '' )157 sn_tmp_rnf = FLD_N( 'runoffs', 24. , 'rotemper' , .TRUE. , .true. , 'yearly' , '' , '' ) 160 158 sn_dep_rnf = FLD_N( 'runoffs', 0. , 'rodepth' , .FALSE. , .true. , 'yearly' , '' , '' ) 161 159 ! … … 218 216 IF ( ln_rnf_att ) THEN 219 217 CALL fld_fill (sf_sal_rnf, (/ sn_sal_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' ) 220 CALL fld_fill (sf_tem_rnf, (/ sn_t em_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )218 CALL fld_fill (sf_tem_rnf, (/ sn_tmp_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' ) 221 219 222 220 rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) … … 248 246 ENDIF 249 247 ! 250 DO jj = 1, jpj251 DO ji = 1, jpi252 rnf_dep(ji,jj) = 0.253 DO jk = 1, rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth254 rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box255 ENDDO256 ENDDO257 ENDDO258 !259 248 260 249 ! ! ======================== -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbc.F90
r2104 r2148 36 36 REAL(wp) :: rn_geoflx_cst = 86.4e-3 ! Constant value of geothermal heat flux 37 37 38 INTEGER , DIMENSION(jpi,jpj) :: nbotlevt ! ocean bottom level index at T-pt39 REAL(wp), DIMENSION(jpi,jpj) :: qgh_trd0 ! geothermal heating trend38 INTEGER , DIMENSION(jpi,jpj) :: nbotlevt ! ocean bottom level index at T-pt 39 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qgh_trd0 ! geothermal heating trend 40 40 41 41 !! * Substitutions -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbl.F90
r2123 r2148 138 138 139 139 140 SUBROUTINE tra_bbl_dif( pt rab, ptraa, kjpt )140 SUBROUTINE tra_bbl_dif( ptb, pta, kjpt ) 141 141 !!---------------------------------------------------------------------- 142 142 !! *** ROUTINE tra_bbl_dif *** … … 155 155 !! convection is satified) 156 156 !! 157 !! ** Action : pt raa increased by the bbl diffusive trend157 !! ** Action : pta increased by the bbl diffusive trend 158 158 !! 159 159 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. … … 161 161 !!---------------------------------------------------------------------- 162 162 INTEGER , INTENT(in ) :: kjpt ! number of tracers 163 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt rab ! before and now tracer fields164 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt raa ! tracer trend163 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 164 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 165 165 !! 166 166 INTEGER :: ji, jj, jn ! dummy loop indices 167 167 INTEGER :: ik ! local integers 168 168 REAL(wp) :: zbtr ! local scalars 169 REAL(wp), DIMENSION(jpi,jpj) :: zptb ! tracer trend 169 170 !!---------------------------------------------------------------------- 170 171 ! 171 ! ! ===========172 172 DO jn = 1, kjpt ! tracer loop 173 173 ! ! =========== 174 # if defined key_vectopt_loop 175 DO jj = 1, 1 ! vector opt. (forced unrolling) 176 DO ji = 1, jpij 177 #else 178 DO jj = 1, jpj 179 DO ji = 1, jpi 180 #endif 181 ik = mbkt(ji,jj) ! bottom T-level index 182 zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S 183 END DO 184 END DO 185 ! ! Compute the trend 174 186 # if defined key_vectopt_loop 175 187 DO jj = 1, 1 ! vector opt. (forced unrolling) … … 181 193 ik = mbkt(ji,jj) ! bottom T-level index 182 194 zbtr = e1e2t_r(ji,jj) / fse3t(ji,jj,ik) 183 pt raa(ji,jj,ik,jn) = ptraa(ji,jj,ik,jn) &184 & + ( ahu_bbl(ji ,jj) * ( ptrab(ji+1,jj ,ik,jn) - ptrab(ji ,jj ,ik,jn) ) &185 & - ahu_bbl(ji ,jj) * ( ptrab(ji ,jj ,ik,jn) - ptrab(ji-1,jj ,ik,jn) ) &186 & + ahv_bbl(ji ,jj) * ( ptrab(ji ,jj+1,ik,jn) - ptrab(ji ,jj ,ik,jn) ) &187 & - ahv_bbl(ji ,jj) * ( ptrab(ji ,jj ,ik,jn) - ptrab(ji ,jj-1,ik,jn) ) ) * zbtr195 pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) & 196 & + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & 197 & - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) & 198 & + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) & 199 & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) * zbtr 188 200 END DO 189 201 END DO … … 194 206 195 207 196 SUBROUTINE tra_bbl_adv( pt rab, ptraa, kjpt )208 SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) 197 209 !!---------------------------------------------------------------------- 198 210 !! *** ROUTINE trc_bbl *** … … 212 224 !!---------------------------------------------------------------------- 213 225 INTEGER , INTENT(in ) :: kjpt ! number of tracers 214 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt rab ! before and now tracer fields215 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt raa ! tracer trend226 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 227 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 216 228 !! 217 229 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 240 252 ! ! up -slope T-point (shelf bottom point) 241 253 zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus) 242 ztra = zu_bbl * ( pt rab(iid,jj,ikus,jn) - ptrab(iis,jj,ikus,jn) ) * zbtr243 pt raa(iis,jj,ikus,jn) = ptraa(iis,jj,ikus,jn) + ztra254 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 255 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 244 256 ! 245 257 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 246 258 zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk) 247 ztra = zu_bbl * ( pt rab(iid,jj,jk+1,jn) - ptrab(iid,jj,jk,jn) ) * zbtr248 pt raa(iid,jj,jk,jn) = ptraa(iid,jj,jk,jn) + ztra259 ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 260 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 249 261 END DO 250 262 ! 251 263 zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud) 252 ztra = zu_bbl * ( pt rab(iis,jj,ikus,jn) - ptrab(iid,jj,ikud,jn) ) * zbtr253 pt raa(iid,jj,ikud,jn) = ptraa(iid,jj,ikud,jn) + ztra264 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 265 pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 254 266 ENDIF 255 267 ! … … 262 274 ! up -slope T-point (shelf bottom point) 263 275 zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs) 264 ztra = zv_bbl * ( pt rab(ji,ijd,ikvs,jn) - ptrab(ji,ijs,ikvs,jn) ) * zbtr265 pt raa(ji,ijs,ikvs,jn) = ptraa(ji,ijs,ikvs,jn) + ztra276 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 277 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 266 278 ! 267 279 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 268 280 zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk) 269 ztra = zv_bbl * ( pt rab(ji,ijd,jk+1,jn) - ptrab(ji,ijd,jk,jn) ) * zbtr270 pt raa(ji,ijd,jk,jn) = ptraa(ji,ijd,jk,jn) + ztra281 ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 282 pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra 271 283 END DO 272 284 ! ! down-slope T-point (deep bottom point) 273 285 zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd) 274 ztra = zv_bbl * ( pt rab(ji,ijs,ikvs,jn) - ptrab(ji,ijd,ikvd,jn) ) * zbtr275 pt raa(ji,ijd,ikvd,jn) = ptraa(ji,ijd,ikvd,jn) + ztra286 ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 287 pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 276 288 ENDIF 277 289 END DO … … 370 382 #endif 371 383 ik = mbkt(ji,jj) ! bottom T-level index 372 ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) ! bottom before T and S373 zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) 384 ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1) ! bottom before T and S 385 zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 374 386 zdep(ji,jj) = fsdept_0(ji,jj,ik) ! bottom T-level reference depth 375 387 ! … … 575 587 ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:) * vmask(:,:,1) 576 588 589 577 590 IF( cp_cfg == "orca" ) THEN !* ORCA configuration : regional enhancement of ah_bbl 578 591 ! -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90
r2120 r2148 15 15 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 16 16 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option 17 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 17 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA 18 !! - ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 18 19 !!---------------------------------------------------------------------- 19 20 20 21 !!---------------------------------------------------------------------- 21 !! tra_nxt : time stepping on t emperature and salinity22 !! tra_nxt_fix : time stepping on t emperature and salinity: fixed volume case23 !! tra_nxt_vvl : time stepping on t emperature and salinity: variable volume case22 !! tra_nxt : time stepping on tracers 23 !! tra_nxt_fix : time stepping on tracers : fixed volume case 24 !! tra_nxt_vvl : time stepping on tracers : variable volume case 24 25 !!---------------------------------------------------------------------- 25 26 USE oce ! ocean dynamics and tracers variables 26 27 USE dom_oce ! ocean space and time domain variables 28 USE sbc_oce ! surface boundary condition: ocean 27 29 USE zdf_oce ! ??? 28 30 USE domvvl ! variable volume … … 38 40 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 39 41 USE prtctl ! Print control 42 USE traqsr ! penetrative solar radiation (needed for nksr) 40 43 USE traswp ! swap array 41 44 USE agrif_opa_update … … 49 52 PUBLIC tra_nxt_vvl ! to be used in trcnxt 50 53 54 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 51 55 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 52 56 … … 96 100 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 97 101 IF(lwp) WRITE(numout,*) '~~~~~~~' 102 ! 103 rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 98 104 ENDIF 99 105 … … 107 113 #endif 108 114 109 #if defined key_obc 115 #if defined key_obc 110 116 IF( lk_obc ) CALL obc_tra( kt ) ! OBC open boundaries 111 117 #endif … … 114 120 #endif 115 121 #if defined key_agrif 116 CALL Agrif_tra 122 CALL Agrif_tra ! AGRIF zoom boundaries 117 123 #endif 118 124 … … 133 139 134 140 ! Leap-Frog + Asselin filter time stepping 135 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts ) ! variable volume level (vvl)136 ELSE ; CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts ) ! fixed volume level141 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts ) ! variable volume level (vvl) 142 ELSE ; CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts ) ! fixed volume level 137 143 ENDIF 138 144 … … 176 182 !! - swap tracer fields to prepare the next time_step. 177 183 !! This can be summurized for tempearture as: 178 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 179 !! ztm = 0 otherwise 184 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 185 !! ztm = 0 otherwise 186 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 180 187 !! tb = tn + atfp*[ tb - 2 tn + ta ] 181 !! tn = ta 188 !! tn = ta 182 189 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 183 190 !! … … 185 192 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 186 193 !!---------------------------------------------------------------------- 187 INTEGER , INTENT(in ) :: kt 188 INTEGER , INTENT(in ) :: kjpt 189 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields190 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields191 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta 194 INTEGER , INTENT(in ) :: kt ! ocean time-step index 195 INTEGER , INTENT(in ) :: kjpt ! number of tracers 196 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 197 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 198 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 192 199 !! 193 200 INTEGER :: ji, jj, jk, jn ! dummy loop indices 194 REAL(wp) :: zt m, ztf! temporary scalars201 REAL(wp) :: ztd, ztm ! temporary scalars 195 202 !!---------------------------------------------------------------------- 196 203 … … 205 212 ! ! (only swap) 206 213 DO jn = 1, kjpt 207 DO jk = 1, jpkm1 214 DO jk = 1, jpkm1 208 215 ptn(:,:,jk,jn) = pta(:,:,jk,jn) ! ptb <-- ptn 209 216 END DO … … 219 226 DO jj = 1, jpj 220 227 DO ji = 1, jpi 221 ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! mean pt 222 ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! Asselin filter on pt 223 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn 224 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 225 pta(ji,jj,jk,jn) = ztm ! pta <-- mean pt 228 ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ! time laplacian on tracers 229 ztm = ptn(ji,jj,jk,jn) + rbcp * ztd ! used for both Asselin and Brown & Campana filters 230 ! 231 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd ! ptb <-- filtered ptn 232 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 233 pta(ji,jj,jk,jn) = ztm ! pta <-- Brown & Campana average 226 234 END DO 227 235 END DO … … 235 243 DO jj = 1, jpj 236 244 DO ji = 1, jpi 237 ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! Asselin filter on t 238 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn 239 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 245 ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ! time laplacian on tracers 246 ! 247 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd ! ptb <-- filtered ptn 248 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 240 249 END DO 241 250 END DO 242 251 END DO 243 252 END DO 244 !245 253 ENDIF 246 254 ! … … 250 258 251 259 252 SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt )260 SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt ) 253 261 !!---------------------------------------------------------------------- 254 262 !! *** ROUTINE tra_nxt_vvl *** … … 263 271 !! - swap tracer fields to prepare the next time_step. 264 272 !! This can be summurized for tempearture as: 265 !! ztm = ( e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T266 !! /( e3t_a +2*e3t_n +e3t_b )267 !! ztm = 0 otherwise273 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 274 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 275 !! ztm = 0 otherwise 268 276 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 269 277 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) … … 274 282 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 275 283 !!---------------------------------------------------------------------- 276 INTEGER , INTENT(in ) :: kt ! ocean time-step index 277 INTEGER , INTENT(in ) :: kjpt ! number of tracers 278 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 279 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 280 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 284 INTEGER , INTENT(in ) :: kt ! ocean time-step index 285 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 286 INTEGER , INTENT(in ) :: kjpt ! number of tracers 287 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 288 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 289 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 281 290 !! 282 INTEGER :: ji, jj, jk, jn ! dummy loop indices 283 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 284 REAL(wp) :: ze3mr, ze3fr ! - - 285 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 291 INTEGER :: ji, jj, jk, jn ! dummy loop indices 292 REAL(wp) :: ztc_a , ztc_n , ztc_b ! temporary scalar 293 REAL(wp) :: ztc_f , ztc_d , ztc_m ! - - 294 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a ! - - 295 REAL(wp) :: ze3t_f, ze3t_d, ze3t_m ! - - 296 REAL :: zfact1, zfact2 ! - - 286 297 !!---------------------------------------------------------------------- 287 298 … … 293 304 ! 294 305 ! 295 IF( neuler == 0 .AND. kt == nit000 ) THEN 296 ! 306 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 307 ! ! (only swap) 297 308 DO jn = 1, kjpt 298 DO jk = 1, jpkm1 299 ptn(:,:,jk,jn) = pta(:,:,jk,jn) 309 DO jk = 1, jpkm1 310 ptn(:,:,jk,jn) = pta(:,:,jk,jn) ! ptb <-- ptn 300 311 END DO 301 312 END DO 302 313 ! 303 ELSE ! general case (Leapfrog + Asselin filter314 ELSE ! general case (Leapfrog + Asselin filter) 304 315 ! 305 ! 306 IF( ln_dynhpg_imp ) THEN 307 ! 308 DO jn = 1, kjpt 316 ! ! ----------------------- ! 317 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 318 ! ! ----------------------- ! 319 DO jn = 1, kjpt 309 320 DO jk = 1, jpkm1 310 321 DO jj = 1, jpj … … 314 325 ze3t_a = fse3t_a(ji,jj,jk) 315 326 ! ! tracer content at Before, now and after 316 ztcb = ptb(ji,jj,jk,jn) * ze3t_b 317 ztcn = ptn(ji,jj,jk,jn) * ze3t_n 318 ztca = pta(ji,jj,jk,jn) * ze3t_a 319 ! 320 ! ! Asselin filter on thickness and tracer content 321 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 322 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 323 ! 324 ! ! filtered tracer including the correction 325 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 326 ztf = ze3fr * ( ztcn + ztc_f ) 327 ! ! mean thickness and tracer 328 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 329 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) 330 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 331 !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr 327 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 328 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 329 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 330 ! ! Time laplacian on tracer contents 331 ! ! used for both Asselin and Brown & Campana filters 332 ze3t_d = ze3t_b - 2.* ze3t_n + ze3t_a 333 ztc_d = ztc_b - 2.* ztc_n + ztc_a 334 ! ! Asselin Filter on thicknesses and tracer contents 335 ztc_f = ztc_n + atfp * ztc_d 336 ztc_m = ztc_n + rbcp * ztc_d 337 ! 338 ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d ) 339 ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d ) 332 340 ! ! swap of arrays 333 ptb(ji,jj,jk,jn) = zt f ! ptb <-- ptn + filter334 pt n(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta335 pt a(ji,jj,jk,jn) = ztm ! pta <-- mean t341 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 342 pta(ji,jj,jk,jn) = ztc_m * ze3t_m ! pta <-- Brown & Campana average 343 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 336 344 END DO 337 345 END DO 338 346 END DO 339 347 END DO 340 ! ! ----------------------- ! 341 ELSE ! explicit hpg case ! 342 ! ! ----------------------- ! 343 DO jn = 1, kjpt 344 DO jk = 1, jpkm1 345 DO jj = 1, jpj 346 DO ji = 1, jpi 347 ze3t_b = fse3t_b(ji,jj,jk) 348 ze3t_n = fse3t_n(ji,jj,jk) 349 ze3t_a = fse3t_a(ji,jj,jk) 350 ! ! tracer content at Before, now and after 351 ztcb = ptb(ji,jj,jk,jn) * ze3t_b 352 ztcn = ptn(ji,jj,jk,jn) * ze3t_n 353 ztca = pta(ji,jj,jk,jn) * ze3t_a 354 ! 355 ! ! Asselin filter on thickness and tracer content 356 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 357 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 358 ! 359 ! ! filtered tracer including the correction 360 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 361 ztf = ( ztcn + ztc_f ) * ze3fr 362 ! ! swap of arrays 363 ptb(ji,jj,jk,jn) = ztf ! tb <-- tn filtered 364 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 348 ! ! ----------------------- ! 349 ELSE ! explicit hpg case ! 350 ! ! ----------------------- ! 351 IF( cdtype == 'TRA' ) THEN 352 ! 353 DO jn = 1, kjpt 354 DO jk = 1, jpkm1 355 zfact1 = atfp * rdttra(jk) 356 zfact2 = zfact1 / rau0 357 DO jj = 1, jpj 358 DO ji = 1, jpi 359 ze3t_b = fse3t_b(ji,jj,jk) 360 ze3t_n = fse3t_n(ji,jj,jk) 361 ze3t_a = fse3t_a(ji,jj,jk) 362 ! ! tracer content at Before, now and after 363 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 364 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 365 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 366 ! 367 ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 368 ztc_f = ztc_n + atfp * ( ztc_a - 2. * ztc_n + ztc_b ) 369 370 IF( jk == 1 ) THEN 371 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 372 IF( jn == jp_tem ) ztc_f = ztc_f - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) ) 373 IF( jn == jp_sal ) ztc_f = ztc_f - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) ) 374 ENDIF 375 IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr ) & 376 & ztc_f = ztc_f - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 377 378 ze3t_f = 1.e0 / ze3t_f 379 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! tb <-- tn filtered 380 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 381 END DO 365 382 END DO 366 383 END DO 367 384 END DO 368 END DO 385 ! 386 ELSE IF( cdtype == 'TRC' ) THEN 387 ! 388 DO jn = 1, kjpt 389 DO jk = 1, jpkm1 390 DO jj = 1, jpj 391 DO ji = 1, jpi 392 ze3t_b = fse3t_b(ji,jj,jk) 393 ze3t_n = fse3t_n(ji,jj,jk) 394 ze3t_a = fse3t_a(ji,jj,jk) 395 ! ! tracer content at Before, now and after 396 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 397 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 398 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 399 ! 400 ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 401 ztc_f = ztc_n + atfp * ( ztc_a - 2. * ztc_n + ztc_b ) 402 403 ze3t_f = 1.e0 / ze3t_f 404 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! tb <-- tn filtered 405 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 406 END DO 407 END DO 408 END DO 409 END DO 410 ! 411 END IF 369 412 ! 370 413 ENDIF -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/traqsr.F90
r2024 r2148 27 27 USE iom ! I/O manager 28 28 USE fldread ! read input fields 29 USE restart ! ocean restart 29 30 30 31 IMPLICIT NONE … … 47 48 ! Module variables 48 49 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 49 INTEGER :: nksr! levels below which the light cannot penetrate ( depth larger than 391 m)50 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m) 50 51 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption 51 52 … … 95 96 REAL(wp) :: zchl, zcoef, zsi0r ! temporary scalars 96 97 REAL(wp) :: zc0, zc1, zc2, zc3 ! - - 98 REAL(wp) :: z1_e3t, zfact ! - - 97 99 REAL(wp), DIMENSION(jpi,jpj) :: zekb, zekg, zekr ! 2D workspace 98 100 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0, ze1 , ze2, ze3, zea ! 3D workspace … … 111 113 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = 0. 112 114 ENDIF 115 116 ! Set before qsr tracer content field 117 ! *********************************** 118 IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1 119 ! ! ----------------------------------- 120 IF( ln_rstart .AND. & ! Restart: read in restart file 121 & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN 122 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field red in the restart file' 123 zfact = 0.5e0 124 CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux 125 ELSE ! No restart or restart not found: Euler forward time stepping 126 zfact = 1.e0 127 qsr_hc_b(:,:,:) = 0.e0 128 ENDIF 129 ELSE ! Swap of forcing field 130 ! ! --------------------- 131 zfact = 0.5e0 132 qsr_hc_b(:,:,:) = qsr_hc_n(:,:,:) 133 ENDIF 134 ! Compute now qsr tracer content field 135 ! ************************************ 113 136 114 137 ! ! ============================================== ! … … 118 141 DO jj = 2, jpjm1 119 142 DO ji = fs_2, fs_jpim1 ! vector opt. 120 ta(ji,jj,jk) = ta(ji,jj,jk) +ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)143 qsr_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 121 144 END DO 122 145 END DO … … 175 198 ! 176 199 DO jk = 1, nksr ! compute and add qsr trend to ta 177 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)200 qsr_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 178 201 END DO 179 202 zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero … … 182 205 ELSE !* Constant Chlorophyll 183 206 DO jk = 1, nksr 184 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) +etot3(:,:,jk) * qsr(:,:)207 qsr_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 185 208 END DO 186 209 ENDIF … … 194 217 DO jj = 2, jpjm1 195 218 DO ji = fs_2, fs_jpim1 ! vector opt. 196 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) +etot3(ji,jj,jk) * qsr(ji,jj)219 qsr_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) 197 220 END DO 198 221 END DO … … 200 223 ! 201 224 ENDIF 225 ! 226 ENDIF 227 ! Add to the general trend 228 ! ************************ 229 DO jk = 1, nksr 230 DO jj = 2, jpjm1 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 z1_e3t = zfact / fse3t(ji,jj,jk) 233 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc_n(ji,jj,jk) ) * z1_e3t 234 END DO 235 END DO 236 END DO 237 ! 238 IF( lrst_oce ) THEN ! Write in the ocean restart file 239 ! ******************************* 240 IF(lwp) WRITE(numout,*) 241 IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ', & 242 & 'at it= ', kt,' date= ', ndastp 243 IF(lwp) WRITE(numout,*) '~~~~' 244 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc_n ) 202 245 ! 203 246 ENDIF -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trasbc.F90
r2113 r2148 7 7 !! 8.2 ! 2001-02 (D. Ludicone) sea ice and free surface 8 8 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! - ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 10 11 !!---------------------------------------------------------------------- 11 12 … … 22 23 USE in_out_manager ! I/O manager 23 24 USE prtctl ! Print control 25 USE restart ! ocean restart 24 26 USE sbcrnf ! River runoff 25 27 USE sbcmod ! ln_rnf 28 USE iom 26 29 27 30 IMPLICIT NONE … … 103 106 INTEGER, INTENT(in) :: kt ! ocean time-step index 104 107 !! 105 INTEGER :: ji, jj, jk ! dummy loop indices 106 REAL(wp) :: zta, zsa ! local scalars, adjustment to temperature and salinity 107 REAL(wp) :: zata, zasa ! local scalars, calculations of automatic change to temp & sal due to vvl (done elsewhere) 108 INTEGER :: ji, jj, jk ! dummy loop indices 109 REAL(wp) :: zta, zsa, zrnf ! local scalars, adjustment to temperature and salinity 108 110 REAL(wp) :: zsrau, zse3t, zdep ! local scalars, 1/density, 1/height of box, 1/height of effected water column 109 111 REAL(wp) :: zdheat, zdsalt ! total change of temperature and salinity 112 REAL(wp) :: zfact, z1_e3t ! 110 113 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 111 114 !!---------------------------------------------------------------------- … … 118 121 119 122 zsrau = 1. / rau0 ! initialization 120 #if defined key_zco121 zse3t = 1. / e3t_0(1)122 #endif123 123 124 124 IF( l_trdtra ) THEN !* Save ta and sa trends … … 127 127 ENDIF 128 128 129 IF( .NOT.ln_traqsr ) qsr(:,:) = 0.e0 ! no solar radiation penetration 130 131 ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff 129 !!gm IF( .NOT.ln_traqsr ) qsr(:,:) = 0.e0 ! no solar radiation penetration 130 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration 131 qns(:,:) = qns(:,:) + qsr(:,:) ! total heat flux in qns 132 qsr(:,:) = 0.e0 ! qsr set to zero 133 ENDIF 134 135 ! Set before sbc tracer content fields 136 ! ************************************ 137 IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1 138 ! ! ----------------------------------- 139 IF( ln_rstart .AND. & ! Restart: read in restart file 140 & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 141 IF(lwp) WRITE(numout,*) ' nit000-1 surface tracer content forcing fields red in the restart file' 142 zfact = 0.5e0 143 CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_hc_b ) ! before heat content sbc trend 144 CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_sc_b ) ! before salt content sbc trend 145 ELSE ! No restart or restart not found: Euler forward time stepping 146 zfact = 1.e0 147 sbc_hc_b(:,:) = 0.e0 148 sbc_sc_b(:,:) = 0.e0 149 ENDIF 150 ELSE ! Swap of forcing fields 151 ! ! ---------------------- 152 zfact = 0.5e0 153 sbc_hc_b(:,:) = sbc_hc_n(:,:) 154 sbc_sc_b(:,:) = sbc_sc_n(:,:) 155 ENDIF 156 ! Compute now sbc tracer content fields 157 ! ************************************* 158 159 ! Concentration dilution effect on (t,s) due to 160 ! evaporation, precipitation and qns, but not river runoff 161 162 IF( lk_vvl ) THEN ! Variable Volume case 163 DO jj = 2, jpj 164 DO ji = fs_2, fs_jpim1 ! vector opt. 165 ! temperature : heat flux + cooling/heating effet of EMP flux 166 sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 167 ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration 168 sbc_sc_n(ji,jj) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) 169 END DO 170 END DO 171 ELSE ! Constant Volume case 172 DO jj = 2, jpj 173 DO ji = fs_2, fs_jpim1 ! vector opt. 174 ! temperature : heat flux 175 sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 176 ! salinity : salt flux + concent./dilut. effect (both in emps) 177 sbc_sc_n(ji,jj) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal) 178 END DO 179 END DO 180 ENDIF 181 ! Concentration dilution effect on (t,s) due to 182 ! river runoff without T, S and depth attributes 183 IF( ln_rnf ) THEN 184 ! 185 IF( lk_vvl ) THEN ! Variable Volume case 186 ! ! cooling/heating effect of runoff & No salinity concent./dilut. effect 187 DO jj = 2, jpj 188 DO ji = fs_2, fs_jpim1 ! vector opt. 189 sbc_hc_n(ji,jj) = sbc_hc_n(ji,jj) + zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_tem) 190 END DO 191 END DO 192 ELSE ! Constant Volume case 193 ! ! concent./dilut. effect only 194 DO jj = 2, jpj 195 DO ji = fs_2, fs_jpim1 ! vector opt. 196 sbc_sc_n(ji,jj) = sbc_sc_n(ji,jj) - zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_sal) 197 END DO 198 END DO 199 ENDIF 200 ! 201 ENDIF 202 ! Add to the general trend 203 ! ************************ 132 204 DO jj = 2, jpj 133 205 DO ji = fs_2, fs_jpim1 ! vector opt. 134 #if ! defined key_zco 135 zse3t = 1. / fse3t(ji,jj,1) 136 #endif 137 IF( lk_vvl) THEN 138 ! temperature : heat flux and heat content of EMP flux 139 zta = ( ro0cpr * qns(ji,jj) - emp(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) ) * zse3t 140 ! Salinity : concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration 141 zsa = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t 142 ELSE 143 zta = ro0cpr * qns(ji,jj) * zse3t ! temperature : heat flux 144 zsa = emps(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t ! salinity : concent./dilut. effect 145 ENDIF 146 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta ! add the trend to the general tracer trend 147 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 206 z1_e3t = zfact / fse3t(ji,jj,1) 207 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + ( sbc_hc_b(ji,jj) + sbc_hc_n(ji,jj) ) * z1_e3t 208 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + ( sbc_sc_b(ji,jj) + sbc_sc_n(ji,jj) ) * z1_e3t 148 209 END DO 149 210 END DO 211 ! Write in the ocean restart file 212 ! ******************************* 213 IF( lrst_oce ) THEN 214 IF(lwp) WRITE(numout,*) 215 IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ', & 216 & 'at it= ', kt,' date= ', ndastp 217 IF(lwp) WRITE(numout,*) '~~~~' 218 CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n ) 219 CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n ) 220 ENDIF 150 221 151 222 IF( ln_rnf .AND. ln_rnf_att ) THEN ! Concentration / dilution effect on (t,s) due to river runoff 223 ! 152 224 DO jj = 1, jpj 153 225 DO ji = 1, jpi 154 zdep = 1. / rnf_dep(ji,jj) 155 zse3t= 1. / fse3t(ji,jj,1) 156 IF( rnf_tem(ji,jj) == -999 ) rnf_tem(ji,jj) = tsn(ji,jj,1,jp_tem) ! if not specified set runoff temp to be sst 226 zdep = 1. / rnf_dep(ji,jj) 227 zse3t = 1. / fse3t(ji,jj,1) 228 rnf_dep(ji,jj) = 0.e0 229 DO jk = 1, rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth 230 rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box 231 END DO 232 IF( rnf_tmp(ji,jj) == -999 ) rnf_tmp(ji,jj) = tsn(ji,jj,1,jp_tem) ! if not specified set runoff temp to be sst 157 233 158 234 IF( rnf(ji,jj) > 0.e0 ) THEN 159 235 ! 236 zrnf = rnf(ji,jj) * zsrau * zdep 160 237 IF( lk_vvl ) THEN 161 238 ! indirect flux, concentration or dilution effect : force a dilution effect in all levels … … 163 240 zdsalt = 0.e0 164 241 DO jk = 1, rnf_mod_dep(ji,jj) 165 zta = -tsn(ji,jj,jk,jp_tem) * rnf(ji,jj) * zsrau * zdep166 zsa = -tsn(ji,jj,jk,jp_sal) * rnf(ji,jj) * zsrau * zdep242 zta = -tsn(ji,jj,jk,jp_tem) * zrnf 243 zsa = -tsn(ji,jj,jk,jp_sal) * zrnf 167 244 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta ! add the trend to the general tracer trend 168 245 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa … … 171 248 END DO 172 249 ! negate this total change in heat and salt content from top level !!gm I don't understand this 173 zta = -zdheat * zse3t 174 zsa = -zdsalt * zse3t 175 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta ! add the trend to the general tracer trend 176 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 177 178 ! direct flux 179 zta = rnf_tem(ji,jj) * rnf(ji,jj) * zsrau * zdep 180 zsa = rnf_sal(ji,jj) * rnf(ji,jj) * zsrau * zdep 250 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) - zdheat * zse3t 251 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) - zdsalt * zse3t 181 252 182 253 DO jk = 1, rnf_mod_dep(ji,jj) 183 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta ! add the trend to the general tracer trend184 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa254 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + rnf_tmp(ji,jj) * zrnf 255 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + rnf_sal(ji,jj) * zrnf 185 256 END DO 186 257 ELSE 187 258 DO jk = 1, rnf_mod_dep(ji,jj) 188 zta = ( rnf_tem(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * rnf(ji,jj) * zsrau * zdep 189 zsa = ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * rnf(ji,jj) * zsrau * zdep 190 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta ! add the trend to the general tracer trend 191 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 259 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( rnf_tmp(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * zrnf 260 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * zrnf 192 261 END DO 193 262 ENDIF 194 263 195 ELSEIF( rnf(ji,jj) < 0.e0 ) THEN ! for use in baltic when flow is out of domain, want no change in temp and sal264 ELSEIF( rnf(ji,jj) < 0.e0 ) THEN ! for use in baltic when flow is out of domain, want no change in temp and sal 196 265 197 266 IF( lk_vvl ) THEN 198 267 ! calculate automatic adjustment to sal and temp due to dilution/concentraion effect 199 zata = tsn(ji,jj,1,jp_tem) * rnf(ji,jj) * zsrau * zse3t 200 zasa = tsn(ji,jj,1,jp_sal) * rnf(ji,jj) * zsrau * zse3t 201 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zata ! add the trend to the general tracer trend 202 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zasa 268 zrnf = rnf(ji,jj) * zsrau * zse3t 269 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * zrnf 270 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * zrnf 203 271 ENDIF 204 272 … … 207 275 END DO 208 276 END DO 209 210 ELSE IF( ln_rnf ) THEN ! Concentration dilution effect on (t,s) due to runoff without T, S and depth attributes 211 212 213 DO jj = 2, jpj 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 #if ! defined key_zco 216 zse3t = 1. / fse3t(ji,jj,1) 217 #endif 218 IF( lk_vvl) THEN 219 zta = rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) * zse3t ! & cooling/heating effect of runoff 220 zsa = 0.e0 ! No salinity concent./dilut. effect 221 ELSE 222 zta = 0.0 ! temperature : heat flux 223 zsa = - rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t ! salinity : concent./dilut. effect 224 ENDIF 225 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta ! add the trend to the general tracer trend 226 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 227 END DO 228 END DO 229 277 ! 230 278 ENDIF 231 279 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/istate.F90
r2104 r2148 67 67 !! ** Purpose : Initialization of the dynamics and tracer fields. 68 68 !!---------------------------------------------------------------------- 69 ! - ML - needed for initialization of e3t_b 70 INTEGER :: jk ! dummy loop indice 69 71 70 72 IF(lwp) WRITE(numout,*) … … 128 130 IF( ln_zps .AND. .NOT. lk_c1d ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! zps: before hor. gradient 129 131 & rhd, gru , grv ) ! of t,s,rd at ocean bottom 130 ! 132 ! 133 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 134 IF( lk_vvl ) THEN 135 DO jk = 1, jpk 136 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 137 ENDDO 138 ENDIF 139 ! 131 140 ENDIF 132 141 ! -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/oce.F90
r2104 r2148 36 36 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 37 37 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m] 38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshf_b , sshf_n , sshf_a!: sea surface height at f-point [m]38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshf_n !: sea surface height at f-point [m] 39 39 ! 40 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: spgu, spgv !: horizontal surface pressure gradient -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/opa.F90
r2104 r2148 281 281 IF( lk_diaar5 ) CALL dia_ar5_init ! ar5 diag 282 282 CALL dia_ptr_init ! Poleward TRansports initialization 283 CALL dia_hsb_init ! heat content, salt content and volume budgets 283 284 CALL trd_mod_init ! Mixed-layer/Vorticity/Integral constraints trends 284 285 ! -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/step.F90
r2104 r2148 221 221 CALL ssh_nxt( kstp ) ! sea surface height at next time step 222 222 223 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 224 223 225 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 224 226 ! Control and restarts -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/step_oce.F90
r2104 r2148 85 85 USE diahth ! thermocline depth (dia_hth routine) 86 86 USE diafwb ! freshwater budget (dia_fwb routine) 87 USE diahsb ! heat, salt and volume budgets (dia_hsb routine) 87 88 USE flo_oce ! floats variables 88 89 USE floats ! floats computation (flo_stp routine) -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/PISCES/p4zmort.F90
r2104 r2148 235 235 !! 236 236 !! ** Method : Read the nampismort namelist and check the parameters 237 !! called at the first timestep (nittrc000)237 !! called at the first timestep 238 238 !! 239 239 !! ** input : Namelist nampismort -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/PISCES/p4zrem.F90
r2104 r2148 408 408 !! 409 409 !! ** Method : Read the nampisrem namelist and check the parameters 410 !! called at the first timestep (nittrc000)410 !! called at the first timestep 411 411 !! 412 412 !! ** input : Namelist nampisrem -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/PISCES/p4zsink.F90
r2104 r2148 321 321 !! 322 322 !! ** Method : Read the nampiskrs namelist and check the parameters 323 !! called at the first timestep (nittrc000)323 !! called at the first timestep 324 324 !! 325 325 !! ** input : Namelist nampiskrs -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcadv.F90
r2082 r2148 65 65 !!---------------------------------------------------------------------- 66 66 67 IF( kt == nit trc000 ) CALL trc_adv_ctl ! initialisation & control of options67 IF( kt == nit000 ) CALL trc_adv_ctl ! initialisation & control of options 68 68 69 IF( neuler == 0 .AND. kt == nit trc000 ) THEN ! at nit00069 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 70 70 r2dt(:) = rdttra(:) * FLOAT(nn_dttrc) ! = rdtra (restarting with Euler time stepping) 71 ELSEIF( kt <= nit trc000 + nn_dttrc ) THEN ! at nit000 or nit000+171 ELSEIF( kt <= nit000 + nn_dttrc ) THEN ! at nit000 or nit000+1 72 72 r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc) ! = 2 rdttra (leapfrog) 73 73 ENDIF -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcdmp.F90
r2030 r2148 82 82 ! 0. Initialization (first time-step only) 83 83 ! -------------- 84 IF( kt == nit trc000 ) CALL trc_dmp_init84 IF( kt == nit000 ) CALL trc_dmp_init 85 85 86 86 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) ! temporary save of trends -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcldf.F90
r2082 r2148 62 62 !!---------------------------------------------------------------------- 63 63 64 IF( kt == nit trc000 ) CALL ldf_ctl ! initialisation & control of options64 IF( kt == nit000 ) CALL ldf_ctl ! initialisation & control of options 65 65 66 66 IF( l_trdtrc ) THEN -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcnxt.F90
r2082 r2148 86 86 !!---------------------------------------------------------------------- 87 87 88 IF( kt == nit trc000 .AND. lwp ) THEN88 IF( kt == nit000 .AND. lwp ) THEN 89 89 WRITE(numout,*) 90 90 WRITE(numout,*) 'trc_nxt : time stepping on passive tracers' … … 109 109 110 110 ! set time step size (Euler/Leapfrog) 111 IF( neuler == 0 .AND. kt == nit trc000) THEN ; r2dt(:) = rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 (Euler)112 ELSEIF( kt <= nit trc000 + 1 ) THEN ; r2dt(:) = 2.* rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 or nit000+1 (Leapfrog)111 IF( neuler == 0 .AND. kt == nit000) THEN ; r2dt(:) = rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 (Euler) 112 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt(:) = 2.* rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 or nit000+1 (Leapfrog) 113 113 ENDIF 114 114 … … 120 120 121 121 ! Leap-Frog + Asselin filter time stepping 122 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nittrc000, trb, trn, tra, jptra ) ! variable volume level (vvl)123 ELSE ; CALL tra_nxt_fix( kt, nittrc000,trb, trn, tra, jptra ) ! fixed volume level122 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, 'TRC', trb, trn, tra, jptra ) ! variable volume level (vvl) 123 ELSE ; CALL tra_nxt_fix( kt, trb, trn, tra, jptra ) ! fixed volume level 124 124 ENDIF 125 125 -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcrad.F90
r2030 r2148 54 54 !!---------------------------------------------------------------------- 55 55 56 IF( kt == nit trc000 ) THEN56 IF( kt == nit000 ) THEN 57 57 IF(lwp) WRITE(numout,*) 58 58 IF(lwp) WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations ' -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcsbc.F90
r2052 r2148 65 65 !! * Local declarations 66 66 INTEGER :: ji, jj, jn ! dummy loop indices 67 REAL(wp) :: z tra, zsrau, zse3t ! temporary scalars67 REAL(wp) :: zsrau, zse3t ! temporary scalars 68 68 REAL(wp), DIMENSION(jpi,jpj) :: zemps ! surface freshwater flux 69 69 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrtrd … … 71 71 !!---------------------------------------------------------------------- 72 72 73 IF( kt == nit trc000 ) THEN73 IF( kt == nit000 ) THEN 74 74 IF(lwp) WRITE(numout,*) 75 75 IF(lwp) WRITE(numout,*) 'trc_sbc : Passive tracers surface boundary condition' … … 82 82 #if ! defined key_offline 83 83 ! Concentration dilution effect on tracer due to evaporation, precipitation, and river runoff 84 IF( lk_vvl ) THEN ; zemps(:,:) = emps(:,:) - emp(:,:) ! volume variable 85 ELSE ; zemps(:,:) = emps(:,:) - rnf(:,:) ! linear free surface 84 IF( lk_vvl ) THEN ! volume variable 85 zemps(:,:) = emps(:,:) - emp(:,:) 86 !!ch zemps(:,:) = 0. 87 ELSE ! linear free surface 88 IF( ln_rnf ) THEN ; zemps(:,:) = emps(:,:) - rnf(:,:) ! E-P-R 89 ELSE ; zemps(:,:) = emps(:,:) 90 ENDIF 86 91 ENDIF 87 92 #else … … 91 96 ENDIF 92 97 #endif 98 93 99 ! 0. initialization 94 100 zsrau = 1. / rau0 95 #if defined key_zco96 zse3t = 1. / e3t_0(1)97 #endif98 99 101 DO jn = 1, jptra 100 102 ! 101 103 IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn) ! save trends 102 104 ! ! add the trend to the general tracer trend 103 105 DO jj = 2, jpj 104 106 DO ji = fs_2, fs_jpim1 ! vector opt. 105 #if ! defined key_zco106 107 zse3t = 1. / fse3t(ji,jj,1) 107 #endif 108 ! add the trend to the general tracer trend 109 ztra = zemps(ji,jj) * zsrau * trn(ji,jj,1,jn) * zse3t 110 tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + ztra 108 tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zemps(ji,jj) * zsrau * trn(ji,jj,1,jn) * zse3t 111 109 END DO 112 110 END DO -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trczdf.F90
r2082 r2148 61 61 !!--------------------------------------------------------------------- 62 62 63 IF( kt == nit trc000 ) CALL zdf_ctl ! initialisation & control of options63 IF( kt == nit000 ) CALL zdf_ctl ! initialisation & control of options 64 64 65 65 #if ! defined key_pisces 66 IF( neuler == 0 .AND. kt == nit trc000 ) THEN ! at nit00066 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 67 67 r2dt(:) = rdttra(:) * FLOAT(nn_dttrc) ! = rdtra (restarting with Euler time stepping) 68 ELSEIF( kt <= nit trc000 + nn_dttrc ) THEN ! at nit000 or nit000+168 ELSEIF( kt <= nit000 + nn_dttrc ) THEN ! at nit000 or nit000+1 69 69 r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc) ! = 2 rdttra (leapfrog) 70 70 ENDIF -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trdmld_trc.F90
r2052 r2148 361 361 !!---------------------------------------------------------------------- 362 362 363 IF( llwarn ) THEN ! warnings 364 IF( ( nittrc000 /= nit000 ) & 365 .OR.( nn_dttrc /= 1 ) ) THEN 366 367 WRITE(numout,*) 'Be careful, trends diags never validated' 368 STOP 'Uncomment this line to proceed' 369 ENDIF 370 ENDIF 363 IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " ) 371 364 372 365 ! ====================================================================== … … 422 415 ! II.1 Set before values of vertically averages passive tracers 423 416 ! ------------------------------------------------------------- 424 IF( kt > nit trc000 ) THEN417 IF( kt > nit000 ) THEN 425 418 DO jn = 1, jptra 426 419 IF( ln_trdtrc(jn) ) THEN … … 507 500 tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc 508 501 509 itmod = kt - nit trc000 + 1502 itmod = kt - nit000 + 1 510 503 it = kt 511 504 … … 915 908 !!---------------------------------------------------------------------- 916 909 ! ... Warnings 917 IF( llwarn ) THEN 918 IF( ( nittrc000 /= nit000 ) & 919 .OR.( nn_dttrc /= 1 ) ) THEN 920 921 WRITE(numout,*) 'Be careful, trends diags never validated' 922 STOP 'Uncomment this line to proceed' 923 END IF 924 END IF 910 IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " ) 925 911 926 912 ! ====================================================================== … … 1038 1024 1039 1025 ! define time axis 1040 itmod = kt - nit trc000 + 11026 itmod = kt - nit000 + 1 1041 1027 it = kt 1042 1028 … … 1318 1304 CALL dia_nam( clhstnam, nn_trd_trc, csuff ) 1319 1305 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & 1320 & 1, jpi, 1, jpj, nit trc000-nn_dttrc, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )1306 & 1, jpi, 1, jpj, nit000, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom ) 1321 1307 1322 1308 !-- Define the ML depth variable … … 1331 1317 CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' ) 1332 1318 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & 1333 & 1, jpi, 1, jpj, nit trc000-nn_dttrc, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )1319 & 1, jpi, 1, jpj, nit000, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom ) 1334 1320 #endif 1335 1321 -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trdmod_trc.F90
r2030 r2148 50 50 !!---------------------------------------------------------------------- 51 51 52 IF( kt == nit trc000 ) THEN52 IF( kt == nit000 ) THEN 53 53 ! IF(lwp)WRITE(numout,*) 54 54 ! IF(lwp)WRITE(numout,*) 'trd_mod_trc:' -
branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/oce_trc.F90
r2104 r2148 201 201 USE sbc_oce , ONLY : emps => emps !: freshwater budget: concentration/dillution [Kg/m2/s] 202 202 USE sbc_oce , ONLY : rnf => rnf !: river runoff [Kg/m2/s] 203 USE sbc_oce , ONLY : ln_rnf => ln_rnf !: runoffs / runoff mouths 203 204 USE sbc_oce , ONLY : fr_i => fr_i !: ice fraction (between 0 to 1) 204 205 USE traqsr , ONLY : rn_abs => rn_abs !: fraction absorbed in the very near surface
Note: See TracChangeset
for help on using the changeset viewer.