Changeset 919 for codes/icosagcm/devel/src/diagnostics
- Timestamp:
- 06/18/19 22:57:38 (5 years ago)
- Location:
- codes/icosagcm/devel/src/diagnostics
- Files:
-
- 1 edited
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/diagnostics/compute_omega.F90
r912 r919 1 MODULE omega_mod 2 1 MODULE compute_omega_mod 3 2 USE icosa 3 IMPLICIT NONE 4 4 PRIVATE 5 5 … … 23 23 END SUBROUTINE W_omega 24 24 25 25 #ifdef BEGIN_DYSL 26 27 KERNEL(compute_omega) 28 29 ! Pressure 30 FORALL_CELLS_EXT() 31 ON_PRIMAL 32 p(CELL) = AP(CELL) + BP(CELL)*ps(HIDX(CELL)) 33 END_BLOCK 34 END_BLOCK 35 36 BARRIER 37 38 ! Mass and grad(ps) 39 FORALL_CELLS_EXT() 40 ON_PRIMAL 41 rhodz(CELL) = (p(CELL)-p(UP(CELL)))*(1./g) 42 END_BLOCK 43 ON_EDGES 44 CST_IFTHEN(IS_BOTTOM_LEVEL) 45 gradps(HIDX(EDGE)) = (ps(HIDX(CELL2))-ps(HIDX(CELL1)))*SIGN*LE 46 CST_ENDIF 47 END_BLOCK 48 END_BLOCK 49 50 ! Mass flux 51 FORALL_CELLS_EXT() 52 ON_EDGES 53 Fe(EDGE)=0.5*(rhodz(CELL1)+rhodz(CELL2)*LE 54 END_BLOCK 55 END_BLOCK 56 57 ! Mass flux divergence 58 ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 59 FORALL_CELLS() 60 ON_PRIMAL 61 divflux=0. 62 FORALL_EDGES 63 divflux = divflux + SIGN*Fe(EDGE) 64 END_BLOCK 65 convm(CELL) = divflux*(1./AI) 66 END_BLOCK 67 END_BLOCK 68 69 ! Barrier needed before and after doing a vertical recurrence 70 BARRIER 71 72 ! vertical integration from up to down 73 SEQUENCE_C1 74 PROLOGUE('llm') 75 convm(CELL)=0. 76 END_BLOCK 77 BODY('llm-1,1,-1') 78 convm(CELL) = convm(CELL) + convm(UP(CELL)) 79 END_BLOCK 80 END_BLOCK 81 82 BARRIER 83 84 ! omega = dp/dt = u.grad p + \pdiff{p}{t} 85 FORALL_CELLS() 86 ON_PRIMAL 87 ugradps=0. 88 FORALL_EDGES 89 ugradps = ugradps + u(EDGE)*gradps(HIDX(EDGE)) 90 END_BLOCK 91 ugradps = .5*(BP(CELL)+BP(UP(CELL)))*ugradps/(-4.*AI) ! sign convention as in Ringler et al. 2010, Eq. 22 p.3072 92 w(CELL) = ugradps - g*.5*(convm(CELL)+convm(UP(CELL))) 93 END_BLOCK 94 END_BLOCK 95 END_BLOCK 96 97 #endif END_DYSL 26 98 27 99 SUBROUTINE compute_omega(ps,u, w) 28 100 USE disvert_mod, ONLY : ap,bp 29 101 USE omp_para 30 IMPLICIT NONE31 102 REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm), ps(iim*jjm) 32 103 REAL(rstd),INTENT(OUT):: w(iim*jjm,llm) 33 104 REAL(rstd):: convm(iim*jjm,llm+1) 34 105 REAL(rstd):: p(iim*jjm,llm+1), rhodz(iim*jjm,llm), Fe(iim*3*jjm,llm) 106 REAL(rstd):: gradps(3*iim*jjm) 35 107 REAL(rstd):: ugradps 36 108 … … 57 129 ENDDO 58 130 ENDDO 131 132 !DIR$ SIMD 133 DO ij=ij_begin_ext, ij_end_ext 134 gradps(ij+u_right) = (ps(ij)-ps(ij+t_right))*ne_right*le(ij+u_right) 135 gradps(ij+u_lup) = (ps(ij)-ps(ij+t_lup)) *ne_lup *le(ij+u_lup) 136 gradps(ij+u_ldown) = (ps(ij)-ps(ij+t_ldown))*ne_ldown*le(ij+u_ldown) 137 END DO 59 138 60 139 !!! Compute mass flux … … 99 178 convm(:,llm+1)=0. 100 179 101 !!! Compute dps102 ! DO j=jj_begin,jj_end103 ! DO i=ii_begin,ii_end104 ! ij=(j-1)*iim+i105 ! ! dps/dt = -int(div flux)dz106 ! dps(ij)=-convm(ij,1) * g107 ! convm(ij,llm+1)=0.108 ! ENDDO109 ! ENDDO110 111 ! Compute Omega = Dp/Dt112 ! with p = A(eta)+B(eta)ps113 ! Dp/Dt = dp/deta.Deta/Dt + B(eta)Dps/Dt114 ! = -mg.Deta/Dt + B.Dps/Dt115 ! By definition the mass flux through model levels is W=m.Deta/Dt with m=-1/g dp/deta116 ! therefore117 ! Dp/Dt = -g.W + B.dps/dt + Bu.grad ps118 ! = B.u.grad ps - g*convm119 120 !!! Compute vertical flux through model layers121 ! DO l = 1,llm-1122 ! DO j=jj_begin,jj_end123 ! DO i=ii_begin,ii_end124 ! ij=(j-1)*iim+i125 ! ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt126 ! ! => w>0 for upward transport127 ! w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) ! g.W = g.convm + B dps/dt128 ! ENDDO129 ! ENDDO130 ! ENDDO131 132 133 180 !!! Compute omega 134 ! -grad ps : ( ne(ij,ldown)*ps(ij,l) + ne(ij+t_ldown,rup)*ps(ij+t_ldown,l) ) ) / de(ij+u_ldown)135 181 136 182 DO l = 1,llm … … 139 185 ij=(j-1)*iim+i 140 186 ugradps = & 141 le(ij+u_right)*u(ij+u_right,l)*( ne(ij,right)*ps(ij) + ne(ij+t_right,left)*ps(ij+t_right) ) & 142 + le(ij+u_rup)*u(ij+u_rup,l)*( ne(ij,rup)*ps(ij) + ne(ij+t_rup,ldown)*ps(ij+t_rup) ) & 143 + le(ij+u_lup)*u(ij+u_lup,l)*( ne(ij,lup)*ps(ij) + ne(ij+t_lup,rdown)*ps(ij+t_lup) ) & 144 + le(ij+u_left)*u(ij+u_left,l)*( ne(ij,left)*ps(ij) + ne(ij+t_left,right)*ps(ij+t_left) ) & 145 + le(ij+u_ldown)*u(ij+u_ldown,l)*( ne(ij,ldown)*ps(ij) + ne(ij+t_ldown,rup)*ps(ij+t_ldown) ) & 146 + le(ij+u_rdown)*u(ij+u_rdown,l)*( ne(ij,rdown)*ps(ij) + ne(ij+t_rdown,lup)*ps(ij+t_rdown) ) 187 u(ij+u_rup,l)*gradps(ij+u_rup) & 188 + u(ij+u_lup,l)*gradps(ij+u_lup) & 189 + u(ij+u_left,l)*gradps(ij+u_left) & 190 + u(ij+u_ldown,l)*gradps(ij+u_ldown) & 191 + u(ij+u_rdown,l)*gradps(ij+u_rdown) & 192 + u(ij+u_right,l)*gradps(ij+u_right) 193 147 194 ugradps = .5*(bp(l)+bp(l+1)) *ugradps/(-4*Ai(ij)) ! sign convention as in Ringler et al. 2010, Eq. 22 p.3072 148 195 w( ij, l) = ugradps - g*.5*(convm( ij,l+1)+convm(ij,l)) … … 155 202 END SUBROUTINE compute_omega 156 203 157 END MODULE omega_mod204 END MODULE compute_omega_mod -
codes/icosagcm/devel/src/diagnostics/observable.f90
r914 r919 48 48 USE vertical_interp_mod 49 49 USE theta2theta_rhodz_mod 50 USE omega_mod50 USE compute_omega_mod, ONLY : w_omega 51 51 USE kinetic_mod 52 52 USE grid_param
Note: See TracChangeset
for help on using the changeset viewer.