![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 138 for codes/icosagcm/trunk/src/advect_tracer.f90
- Timestamp:
- 02/16/13 17:03:57 (11 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/advect_tracer.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect_tracer.f90
r137 r138 1 1 MODULE advect_tracer_mod 2 2 USE icosa 3 IMPLICIT NONE 3 4 PRIVATE 4 5 … … 8 9 TYPE(t_field),POINTER :: f_cc(:) ! starting point of backward-trajectory (Miura approach) 9 10 10 ! TYPE(t_field),POINTER :: f_rhodzm1(:)11 ! TYPE(t_field),POINTER :: f_rhodz(:)12 13 11 REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 14 12 … … 19 17 SUBROUTINE init_advect_tracer 20 18 USE advect_mod 21 IMPLICIT NONE22 19 REAL(rstd),POINTER :: tangent(:,:) 23 20 REAL(rstd),POINTER :: normal(:,:) 24 21 INTEGER :: ind 25 22 26 CALL allocate_field(f_normal,field_u,type_real,3) 27 CALL allocate_field(f_tangent,field_u,type_real,3) 28 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 29 CALL allocate_field(f_cc,field_u,type_real,llm,3) 30 31 ! CALL allocate_field(f_rhodz,field_t,type_real,llm) 32 ! CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 23 CALL allocate_field(f_normal,field_u,type_real,3, name='normal') 24 CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent') 25 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d') 26 CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc') 33 27 34 28 DO ind=1,ndomain … … 43 37 44 38 SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 45 USE icosa46 USE field_mod47 USE domain_mod48 USE dimensions49 USE grid_param50 USE geometry51 USE metric52 39 USE advect_mod 53 USE advect_mod54 USE disvert_mod55 40 USE mpipara 56 41 IMPLICIT NONE … … 62 47 TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step 63 48 64 REAL(rstd),POINTER :: q(:,:,:), normal(:,: ,:), tangent(:,:,:), gradq3d(:,:,:)49 REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), gradq3d(:,:,:), cc(:,:,:) 65 50 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 66 51 REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 67 INTEGER :: ind 52 INTEGER :: ind,k 68 53 69 54 CALL transfert_request(f_u,req_e1) 70 CALL transfert_request(f_u,req_e1) 55 ! CALL transfert_request(f_hfluxt,req_e1) ! BUG : This (unnecessary) transfer makes the computation go wrong 56 CALL transfert_request(f_wfluxt,req_i1) 71 57 CALL transfert_request(f_q,req_i1) 58 CALL transfert_request(f_rhodz,req_i1) 72 59 73 60 IF (is_mpi_root) PRINT *, 'Advection scheme' 74 61 62 ! DO ind=1,ndomain 63 ! CALL swap_dimensions(ind) 64 ! CALL swap_geometry(ind) 65 ! normal = f_normal(ind) 66 ! tangent = f_tangent(ind) 67 ! cc = f_cc(ind) 68 ! u = f_u(ind) 69 ! q = f_q(ind) 70 ! rhodz = f_rhodz(ind) 71 ! hfluxt = f_hfluxt(ind) 72 ! wfluxt = f_wfluxt(ind) 73 ! gradq3d = f_gradq3d(ind) 74 ! 75 ! ! 1/2 vertical transport 76 ! DO k = 1, nqtot 77 ! CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 78 ! END DO 79 ! 80 ! ! horizontal transport 81 ! CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 82 ! DO k = 1,nqtot 83 ! CALL compute_gradq3d(q(:,:,k),gradq3d) 84 ! CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 85 ! END DO 86 ! 87 ! ! 1/2 vertical transport 88 ! DO k = 1,nqtot 89 ! CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 90 ! END DO 91 ! END DO 92 93 ! 1/2 vertical transport + back-trajectories 94 DO ind=1,ndomain 95 CALL swap_dimensions(ind) 96 CALL swap_geometry(ind) 97 normal = f_normal(ind) 98 tangent = f_tangent(ind) 99 cc = f_cc(ind) 100 u = f_u(ind) 101 q = f_q(ind) 102 rhodz = f_rhodz(ind) 103 wfluxt = f_wfluxt(ind) 104 DO k = 1, nqtot 105 CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 106 END DO 107 CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 108 END DO 109 110 CALL transfert_request(f_q,req_i1) ! necessary ? 111 CALL transfert_request(f_rhodz,req_i1) ! necessary ? 112 113 ! horizontal transport - split in two to place transfer of gradq3d 114 115 DO k = 1, nqtot 116 DO ind=1,ndomain 117 CALL swap_dimensions(ind) 118 CALL swap_geometry(ind) 119 q = f_q(ind) 120 gradq3d = f_gradq3d(ind) 121 CALL compute_gradq3d(q(:,:,k),gradq3d) 122 END DO 123 124 CALL transfert_request(f_gradq3d,req_i1) 125 126 DO ind=1,ndomain 127 CALL swap_dimensions(ind) 128 CALL swap_geometry(ind) 129 cc = f_cc(ind) 130 q = f_q(ind) 131 rhodz = f_rhodz(ind) 132 hfluxt = f_hfluxt(ind) 133 gradq3d = f_gradq3d(ind) 134 CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 135 END DO 136 END DO 137 138 CALL transfert_request(f_q,req_i1) ! necessary ? 139 CALL transfert_request(f_rhodz,req_i1) ! necessary ? 140 141 ! 1/2 vertical transport 75 142 DO ind=1,ndomain 76 143 CALL swap_dimensions(ind) … … 78 145 q = f_q(ind) 79 146 rhodz = f_rhodz(ind) 80 hfluxt = f_hfluxt(ind)81 147 wfluxt = f_wfluxt(ind) 82 normal = f_normal(ind) 83 tangent = f_tangent(ind) 84 gradq3d = f_gradq3d(ind) 85 CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 148 DO k = 1,nqtot 149 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 150 END DO 86 151 END DO 87 152 88 153 END SUBROUTINE advect_tracer 89 154 90 SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d)91 USE icosa92 USE field_mod93 USE domain_mod94 USE dimensions95 USE grid_param96 USE geometry97 USE metric98 USE advect_mod99 USE advect_mod100 USE disvert_mod101 USE mpipara102 REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:)103 REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:)104 REAL(rstd), INTENT(INOUT) :: rhodz(:,:)105 REAL(rstd), INTENT(INOUT) :: q(:,:,:)106 REAL(rstd), INTENT(OUT) :: gradq3d(:,:,:)107 INTEGER :: k108 ! for mass/tracer consistency each transport step also updates the mass field rhodz109 ! mass is updated after all tracers have been updated, when k==nqtot110 111 ! 1/2 vertical transport112 DO k = 1, nqtot113 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k))114 END DO115 116 ! horizontal transport117 DO k = 1,nqtot118 CALL compute_gradq3d(q(:,:,k),gradq3d)119 CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv)120 END DO121 122 ! 1/2 vertical transport123 DO k = 1,nqtot124 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))125 END DO126 END SUBROUTINE compute_adv_tracer127 128 155 SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q) 129 156 ! … … 135 162 ! wfluxt >0 for upward transport 136 163 ! ******************************************************************** 137 USE icosa138 164 IMPLICIT NONE 139 165 LOGICAL, INTENT(IN) :: update_mass … … 152 178 ! finite difference of q 153 179 DO l=2,llm 154 DO j=jj_begin ,jj_end155 DO i=ii_begin ,ii_end180 DO j=jj_begin-1,jj_end+1 181 DO i=ii_begin-1,ii_end+1 156 182 ij=(j-1)*iim+i 157 183 dzqw(ij,l)=q(ij,l)-q(ij,l-1) … … 164 190 ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 165 191 DO l=2,llm-1 166 DO j=jj_begin ,jj_end167 DO i=ii_begin ,ii_end192 DO j=jj_begin-1,jj_end+1 193 DO i=ii_begin-1,ii_end+1 168 194 ij=(j-1)*iim+i 169 195 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN … … 179 205 180 206 ! 0 slope in top and bottom layers 181 DO j=jj_begin ,jj_end182 DO i=ii_begin ,ii_end207 DO j=jj_begin-1,jj_end+1 208 DO i=ii_begin-1,ii_end+1 183 209 ij=(j-1)*iim+i 184 210 dzq(ij,1)=0. … … 190 216 ! then amount of q leaving level l/l+1 = wq = w * qq 191 217 DO l = 1,llm-1 192 DO j=jj_begin,jj_end 193 DO i=ii_begin,ii_end 194 ij=(j-1)*iim+i 195 IF(wfluxt(ij,l+1).gt.0.) THEN 196 ! downward transport, sigw<0 197 ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 198 w = fac*wfluxt(ij,l+1) 218 DO j=jj_begin-1,jj_end+1 219 DO i=ii_begin-1,ii_end+1 220 ij=(j-1)*iim+i 221 w = fac*wfluxt(ij,l+1) 222 IF(w>0.) THEN ! upward transport, upwind side is at level l 223 sigw = w/mass(ij,l) 224 qq = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 225 ELSE ! downward transport, upwind side is at level l+1 199 226 sigw = w/mass(ij,l+1) 200 qq = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1) 201 wq(ij,l+1) = w*qq 202 ELSE 203 ! upward transport, sigw>0 204 ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 205 w = fac*wfluxt(ij,l) 206 sigw = w/mass(ij,l) 207 qq = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) 208 wq(ij,l+1) = w*qq 227 qq = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 209 228 ENDIF 229 wq(ij,l+1) = w*qq 210 230 ENDDO 211 231 ENDDO 212 232 END DO 213 233 ! wq = 0 at top and bottom 214 DO j=jj_begin ,jj_end215 DO i=ii_begin ,ii_end234 DO j=jj_begin-1,jj_end+1 235 DO i=ii_begin-1,ii_end+1 216 236 ij=(j-1)*iim+i 217 237 wq(ij,llm+1)=0. … … 222 242 ! update q, mass is updated only after all q's have been updated 223 243 DO l=1,llm 224 DO j=jj_begin ,jj_end225 DO i=ii_begin ,ii_end244 DO j=jj_begin-1,jj_end+1 245 DO i=ii_begin-1,ii_end+1 226 246 ij=(j-1)*iim+i 227 247 newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1))
Note: See TracChangeset
for help on using the changeset viewer.
