Changeset 599 for codes/icosagcm/trunk
- Timestamp:
- 10/19/17 17:04:26 (7 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 5 added
- 1 deleted
- 8 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/diagnostics/check_conserve.f90
r548 r599 299 299 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv 300 300 301 REAL(rstd) :: ucenter(iim*jjm, 3,llm)301 REAL(rstd) :: ucenter(iim*jjm,llm,3) 302 302 REAL(rstd) :: ulon(iim*jjm,llm) 303 303 REAL(rstd) :: ulat(iim*jjm,llm) -
codes/icosagcm/trunk/src/diagnostics/observable.f90
r581 r599 6 6 TYPE(t_field),POINTER, SAVE :: f_buf_i(:), & 7 7 f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH 8 f_buf_ulon(:), f_buf_ulat(:), & 9 f_buf_u3d(:) ! unused, remove ? 8 f_buf_ulon(:), f_buf_ulat(:) 10 9 TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:) 11 10 TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:) … … 24 23 CALL allocate_field(f_buf2_i, field_t,type_real,llm,name="buffer2_i") 25 24 CALL allocate_field(f_buf_p, field_t,type_real,llm+1) 26 CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm) ! 3D vel at cell centers27 25 CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon") 28 26 CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat") … … 49 47 USE theta2theta_rhodz_mod 50 48 USE omega_mod 49 USE diagflux_mod 51 50 LOGICAL, INTENT(IN) :: init 52 51 INTEGER :: l … … 160 159 END IF 161 160 161 IF(.NOT. init) THEN 162 IF(diagflux_on) THEN 163 CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat) 164 CALL output_field("mass_t", f_masst) 165 CALL output_field("massflux_lon",f_buf_ulon) 166 CALL output_field("massflux_lat",f_buf_ulat) 167 168 CALL transfert_request(f_epotfluxt,req_e1_vect) 169 CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat) 170 CALL output_field("epot_t", f_epot) 171 CALL output_field("epotflux_lon",f_buf_ulon) 172 CALL output_field("epotflux_lat",f_buf_ulat) 173 174 CALL transfert_request(f_ekinfluxt,req_e1_vect) 175 CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat) 176 CALL output_field("ekin_t", f_ekin) 177 CALL output_field("ekinflux_lon",f_buf_ulon) 178 CALL output_field("ekinflux_lat",f_buf_ulat) 179 180 CALL transfert_request(f_enthalpyfluxt,req_e1_vect) 181 CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat) 182 CALL output_field("enthalpy_t", f_enthalpy) 183 CALL output_field("enthalpyflux_lon",f_buf_ulon) 184 CALL output_field("enthalpyflux_lat",f_buf_ulat) 185 186 CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat) 187 CALL output_field("qmass_t", f_qmasst) 188 CALL output_field("qflux_lon",f_qfluxt_lon) 189 CALL output_field("qflux_lat",f_qfluxt_lat) 190 CALL zero_qfluxt ! restart accumulating fluxes 191 END IF 192 END IF 162 193 END SUBROUTINE write_output_fields_basic 163 194 -
codes/icosagcm/trunk/src/diagnostics/wind.F90
r581 r599 1 1 MODULE wind_mod 2 3 CONTAINS 4 5 SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat) 2 USE omp_para 6 3 USE icosa 7 4 IMPLICIT NONE 5 PRIVATE 6 7 PUBLIC :: compute_wind_centered, compute_flux_centered, & 8 compute_wind_centered_lonlat_compound, compute_wind2d_perp_from_lonlat_compound, & 9 compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound, & 10 un2ulonlat, ulonlat2un 11 12 CONTAINS 13 14 SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat) 8 15 TYPE(t_field), POINTER :: f_u(:) ! IN : normal velocity components on edges 9 TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons 10 16 TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons 11 17 REAL(rstd),POINTER :: u(:,:), ulon(:,:), ulat(:,:) 12 18 INTEGER :: ind … … 25 31 26 32 SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u) 27 USE icosa28 IMPLICIT NONE29 33 TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons 30 34 TYPE(t_field), POINTER :: f_u(:) ! OUT : normal velocity components on edges … … 46 50 47 51 SUBROUTINE compute_wind_centered(ue,ucenter) 48 USE icosa49 USE omp_para50 IMPLICIT NONE51 52 REAL(rstd) :: ue(3*iim*jjm,llm) 52 REAL(rstd) :: ucenter(iim*jjm,3,llm) 53 INTEGER :: i,j,ij,l 54 53 REAL(rstd) :: ucenter(iim*jjm,llm,3) 54 INTEGER :: ij,l 55 REAL(rstd), PARAMETER :: scale=1. 56 REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz 57 #include "../kernels/wind_centered.k90" 58 END SUBROUTINE compute_wind_centered 59 60 SUBROUTINE compute_flux_centered(scale,ue,ucenter) 61 REAL(rstd), INTENT(IN) :: scale 62 REAL(rstd) :: ue(3*iim*jjm,llm) 63 REAL(rstd) :: ucenter(iim*jjm,llm,3) 64 INTEGER :: ij,l 65 REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz 66 #include "../kernels/flux_centered.k90" 67 END SUBROUTINE compute_flux_centered 68 69 70 SUBROUTINE compute_wind_on_edge(ue,uedge) 71 REAL(rstd) :: ue(3*iim*jjm,llm) 72 REAL(rstd) :: uedge(3*iim*jjm,llm,3) 73 74 REAL(rstd) :: ut(3*iim*jjm,llm) 75 INTEGER :: i,j,ij,l 76 77 CALL compute_tangential_compound(ue,ut) 78 55 79 DO l=ll_begin,ll_end 56 80 DO j=jj_begin,jj_end 57 81 DO i=ii_begin,ii_end 58 82 ij=(j-1)*iim+i 59 ucenter(ij,:,l)=1/Ai(ij)* & 60 ( ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:)) & 61 + ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:)) & 62 + ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:)) & 63 + ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:)) & 64 + ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))& 65 + ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:))) 66 ENDDO 67 ENDDO 68 ENDDO 69 70 END SUBROUTINE compute_wind_centered 71 72 73 SUBROUTINE compute_wind_on_edge(ue,uedge) 74 USE icosa 75 USE omp_para 76 77 IMPLICIT NONE 78 REAL(rstd) :: ue(3*iim*jjm,llm) 79 REAL(rstd) :: uedge(3*iim*jjm,3,llm) 80 81 REAL(rstd) :: ut(3*iim*jjm,llm) 82 INTEGER :: i,j,ij,l 83 84 CALL compute_tangential_compound(ue,ut) 85 86 DO l=ll_begin,ll_end 87 DO j=jj_begin,jj_end 88 DO i=ii_begin,ii_end 89 ij=(j-1)*iim+i 90 uedge(ij+u_right,:,l)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right) 91 uedge(ij+u_lup,:,l)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup) 92 uedge(ij+u_ldown,:,l)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown) 83 uedge(ij+u_right,l,:)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right) 84 uedge(ij+u_lup,l,:)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup) 85 uedge(ij+u_ldown,l,:)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown) 93 86 ENDDO 94 87 ENDDO … … 100 93 101 94 SUBROUTINE compute_tangential_compound(ue,ut) 102 USE icosa103 USE omp_para104 IMPLICIT NONE105 95 REAL(rstd) :: ue(3*iim*jjm,llm) 106 96 REAL(rstd) :: ut(3*iim*jjm,llm) … … 155 145 END SUBROUTINE compute_tangential_compound 156 146 157 SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat) 158 USE icosa 159 USE omp_para 160 161 IMPLICIT NONE 162 REAL(rstd) :: u(3*iim*jjm,3,llm) 163 REAL(rstd) :: ulon(3*iim*jjm,3,llm) 164 REAL(rstd) :: ulat(3*iim*jjm,3,llm) 165 166 INTEGER :: i,j,ij,l 167 147 ! SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat) 148 ! REAL(rstd) :: u(3*iim*jjm,3,llm) 149 ! REAL(rstd) :: ulon(3*iim*jjm,3,llm) 150 ! REAL(rstd) :: ulat(3*iim*jjm,3,llm) 151 ! 152 ! INTEGER :: i,j,ij,l 153 ! 154 ! 155 ! DO l=ll_begin,ll_end 156 ! DO j=jj_begin-1,jj_end+1 157 ! DO i=ii_begin-1,ii_end+1 158 ! ij=(j-1)*iim+i 159 ! ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:) 160 ! ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:) 161 ! ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:) 162 ! 163 ! ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:) 164 ! ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:) 165 ! ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:) 166 ! 167 ! ENDDO 168 ! ENDDO 169 ! ENDDO 170 ! 171 ! END SUBROUTINE compute_wind_lonlat_compound 172 173 SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u) 174 REAL(rstd) :: u(3*iim*jjm,llm,3) 175 REAL(rstd) :: ulon(3*iim*jjm,llm) 176 REAL(rstd) :: ulat(3*iim*jjm,llm) 177 178 INTEGER :: i,j,ij,l 168 179 169 180 DO l=ll_begin,ll_end … … 171 182 DO i=ii_begin-1,ii_end+1 172 183 ij=(j-1)*iim+i 173 ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:) 174 ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:) 175 ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:) 176 177 ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:) 178 ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:) 179 ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:) 180 181 ENDDO 182 ENDDO 183 ENDDO 184 185 END SUBROUTINE compute_wind_lonlat_compound 186 187 SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u) 188 USE icosa 189 USE omp_para 190 191 IMPLICIT NONE 192 REAL(rstd) :: u(3*iim*jjm,3,llm) 193 REAL(rstd) :: ulon(3*iim*jjm,llm) 194 REAL(rstd) :: ulat(3*iim*jjm,llm) 195 196 INTEGER :: i,j,ij,l 197 198 DO l=ll_begin,ll_end 199 DO j=jj_begin-1,jj_end+1 200 DO i=ii_begin-1,ii_end+1 201 ij=(j-1)*iim+i 202 u(ij+u_right,:,l)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:) 203 u(ij+u_lup,:,l)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:) 204 u(ij+u_ldown,:,l)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:) 184 u(ij+u_right,l,:)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:) 185 u(ij+u_lup,l,:)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:) 186 u(ij+u_ldown,l,:)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:) 205 187 ENDDO 206 188 ENDDO … … 210 192 211 193 SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u) 212 USE icosa 213 USE omp_para 214 215 IMPLICIT NONE 216 REAL(rstd) :: u(iim*jjm,3,llm) 194 REAL(rstd) :: u(iim*jjm,llm,3) 217 195 REAL(rstd) :: ulon(iim*jjm,llm) 218 196 REAL(rstd) :: ulat(iim*jjm,llm) 219 220 197 INTEGER :: i,j,ij,l 221 198 DO l=ll_begin,ll_end … … 223 200 DO i=ii_begin-1,ii_end+1 224 201 ij=(j-1)*iim+i 225 u(ij,:,l)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:) 226 ENDDO 227 ENDDO 228 ENDDO 229 202 u(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:) 203 ENDDO 204 ENDDO 205 ENDDO 230 206 END SUBROUTINE compute_wind_centered_from_lonlat_compound 231 207 232 208 SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u) 233 USE icosa234 235 IMPLICIT NONE236 209 REAL(rstd) :: u(3*iim*jjm,3) 237 210 REAL(rstd) :: ulon(3*iim*jjm) … … 252 225 253 226 SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up) 254 USE icosa255 USE omp_para256 257 IMPLICIT NONE258 227 REAL(rstd) :: up(3*iim*jjm,llm) 259 228 REAL(rstd) :: ulon(3*iim*jjm,llm) 260 229 REAL(rstd) :: ulat(3*iim*jjm,llm) 261 REAL(rstd) :: u(3*iim*jjm, 3,llm)230 REAL(rstd) :: u(3*iim*jjm,llm,3) 262 231 263 232 INTEGER :: i,j,ij,l … … 269 238 DO i=ii_begin-1,ii_end+1 270 239 ij=(j-1)*iim+i 271 up(ij+u_right,l)=sum(u(ij+u_right, :,l)*ep_e(ij+u_right,:))272 up(ij+u_lup,l)=sum(u(ij+u_lup, :,l)*ep_e(ij+u_lup,:))273 up(ij+u_ldown,l)=sum(u(ij+u_ldown, :,l)*ep_e(ij+u_ldown,:))240 up(ij+u_right,l)=sum(u(ij+u_right,l,:)*ep_e(ij+u_right,:)) 241 up(ij+u_lup,l)=sum(u(ij+u_lup,l,:)*ep_e(ij+u_lup,:)) 242 up(ij+u_ldown,l)=sum(u(ij+u_ldown,l,:)*ep_e(ij+u_ldown,:)) 274 243 ENDDO 275 244 ENDDO … … 279 248 280 249 SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up) 281 USE icosa282 283 IMPLICIT NONE284 250 REAL(rstd) :: up(3*iim*jjm) 285 251 REAL(rstd) :: ulon(3*iim*jjm) … … 302 268 303 269 SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat) 304 USE icosa 305 USE omp_para 306 307 IMPLICIT NONE 308 REAL(rstd) :: uc(iim*jjm,3,llm) 270 REAL(rstd) :: uc(iim*jjm,llm,3) 309 271 REAL(rstd) :: ulon(iim*jjm,llm) 310 272 REAL(rstd) :: ulat(iim*jjm,llm) 311 273 312 274 INTEGER :: i,j,ij,l 313 314 275 315 276 DO l=ll_begin,ll_end … … 317 278 DO i=ii_begin,ii_end 318 279 ij=(j-1)*iim+i 319 ulon(ij,l)=sum(uc(ij, :,l)*elon_i(ij,:))320 ulat(ij,l)=sum(uc(ij, :,l)*elat_i(ij,:))280 ulon(ij,l)=sum(uc(ij,l,:)*elon_i(ij,:)) 281 ulat(ij,l)=sum(uc(ij,l,:)*elat_i(ij,:)) 321 282 ENDDO 322 283 ENDDO … … 326 287 327 288 SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc) 328 USE icosa329 USE omp_para330 331 IMPLICIT NONE332 289 REAL(rstd) :: ulon(iim*jjm,llm) 333 290 REAL(rstd) :: ulat(iim*jjm,llm) 334 REAL(rstd) :: uc(iim*jjm, 3,llm)291 REAL(rstd) :: uc(iim*jjm,llm,3) 335 292 336 293 INTEGER :: i,j,ij,l … … 341 298 DO i=ii_begin,ii_end 342 299 ij=(j-1)*iim+i 343 uc(ij, :,l)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:)300 uc(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:) 344 301 ENDDO 345 302 ENDDO … … 348 305 END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered 349 306 350 351 352 307 SUBROUTINE compute_wind_perp_from_wind_centered(uc,un) 353 USE icosa 354 USE omp_para 355 308 356 309 IMPLICIT NONE 357 REAL(rstd),INTENT(IN) :: uc(iim*jjm, 3,llm)310 REAL(rstd),INTENT(IN) :: uc(iim*jjm,llm,3) 358 311 REAL(rstd),INTENT(OUT) :: un(3*iim*jjm,llm) 359 312 … … 365 318 DO i=ii_begin,ii_end 366 319 ij=(j-1)*iim+i 367 un(ij+u_right,l) = sum(0.5*(uc(ij, :,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:))368 un(ij+u_lup,l) = sum(0.5*(uc(ij, :,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:))369 un(ij+u_ldown,l) = sum(0.5*(uc(ij, :,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:))320 un(ij+u_right,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:)) 321 un(ij+u_lup,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:)) 322 un(ij+u_ldown,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:)) 370 323 ENDDO 371 324 ENDDO … … 376 329 377 330 SUBROUTINE compute_un2ulonlat(un, ulon, ulat) 378 USE icosa379 380 IMPLICIT NONE381 331 REAL(rstd),INTENT(IN) :: un(3*iim*jjm,llm) 382 332 REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm) 383 333 REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm) 384 334 385 REAL(rstd) :: uc(iim*jjm, 3,llm)335 REAL(rstd) :: uc(iim*jjm,llm,3) 386 336 387 337 CALL compute_wind_centered(un,uc) … … 391 341 392 342 SUBROUTINE compute_ulonlat2un(ulon, ulat,un) 393 USE icosa394 395 IMPLICIT NONE396 343 REAL(rstd),INTENT(IN) :: ulon(iim*jjm,llm) 397 344 REAL(rstd),INTENT(IN) :: ulat(iim*jjm,llm) 398 345 REAL(rstd),INTENT(OUT) :: un(3*iim*jjm,llm) 399 346 400 REAL(rstd) :: uc(iim*jjm, 3,llm)347 REAL(rstd) :: uc(iim*jjm,llm,3) 401 348 402 349 CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc) -
codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90
r580 r599 5 5 USE caldyn_kernels_base_mod 6 6 USE caldyn_kernels_mod 7 USE omp_para 8 USE mpipara 7 9 IMPLICIT NONE 8 10 PRIVATE … … 13 15 14 16 SUBROUTINE init_caldyn 15 USE icosa16 USE observable_mod17 USE mpipara18 USE omp_para19 IMPLICIT NONE20 17 CHARACTER(len=255) :: def 21 18 INTEGER :: ind … … 124 121 125 122 SUBROUTINE allocate_caldyn 126 USE icosa127 IMPLICIT NONE128 129 123 CALL allocate_field(f_out_u,field_u,type_real,llm) 130 124 CALL allocate_field(f_qu,field_u,type_real,llm) … … 142 136 143 137 SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux) 144 USE icosa145 USE mpipara146 USE omp_para147 138 TYPE(t_field),POINTER :: f_phis(:) 148 139 TYPE(t_field),POINTER :: f_geopot(:) … … 186 177 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 187 178 f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 188 USE icosa189 179 USE observable_mod 190 180 USE disvert_mod, ONLY : caldyn_eta, eta_mass 191 USE vorticity_mod192 USE kinetic_mod193 USE theta2theta_rhodz_mod194 USE wind_mod195 USE mpipara196 181 USE trace 197 USE omp_para198 USE output_field_mod199 USE checksum_mod200 IMPLICIT NONE201 182 LOGICAL,INTENT(IN) :: write_out 202 183 TYPE(t_field),POINTER :: f_phis(:) … … 375 356 376 357 SUBROUTINE check_mass_conservation(f_ps,f_dps) 377 USE icosa378 USE mpipara379 IMPLICIT NONE380 358 TYPE(t_field),POINTER :: f_ps(:) 381 359 TYPE(t_field),POINTER :: f_dps(:) -
codes/icosagcm/trunk/src/icosa_init.f90
r554 r599 1 1 MODULE icosa_init_mod 2 3 4 2 5 3 CONTAINS … … 22 20 USE restart_mod 23 21 USE etat0_mod 22 USE diagflux_mod 24 23 IMPLICIT NONE 25 24 … … 55 54 CALL init_timeloop 56 55 CALL init_physics 57 56 57 CALL init_diagflux 58 58 CALL timeloop 59 59 CALL switch_omp_no_distrib_level -
codes/icosagcm/trunk/src/physics/physics.f90
r548 r599 247 247 248 248 REAL(rstd) :: p(iim*jjm,llm+1) 249 REAL(rstd) :: uc(iim*jjm, 3,llm)249 REAL(rstd) :: uc(iim*jjm,llm,3) 250 250 REAL(rstd) :: ulon(iim*jjm,llm) 251 251 REAL(rstd) :: ulat(iim*jjm,llm) … … 311 311 REAL(rstd) :: dulat(iim*jjm,llm) 312 312 REAL(rstd) :: ue(3*iim*jjm,llm) 313 REAL(rstd) :: duc(iim*jjm, 3,llm)313 REAL(rstd) :: duc(iim*jjm,llm,3) 314 314 REAL(rstd) :: dt2, due 315 315 INTEGER :: i,j,ij,l … … 321 321 DO i=ii_begin,ii_end 322 322 ij=(j-1)*iim+i 323 due = sum( (duc(ij, :,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )323 due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) ) 324 324 ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due 325 325 326 due = sum( (duc(ij, :,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )326 due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) ) 327 327 ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due 328 328 329 due = sum( (duc(ij, :,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )329 due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) ) 330 330 ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due 331 331 ENDDO -
codes/icosagcm/trunk/src/physics/physics_lmdz_generic.F90
r548 r599 78 78 CALL allocate_field(f_dq,field_t,type_real,llm,nqtot) 79 79 CALL allocate_field(f_dps,field_t,type_real) 80 CALL allocate_field(f_duc,field_t,type_real, 3,llm)80 CALL allocate_field(f_duc,field_t,type_real,llm,3) 81 81 !$OMP END PARALLEL 82 82 … … 362 362 DO i=ii_begin,ii_end 363 363 ij=(j-1)*iim+i 364 duc(ij, :,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)364 duc(ij,l,:)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:) 365 365 ENDDO 366 366 ENDDO … … 371 371 DO i=ii_begin,ii_end 372 372 ij=(j-1)*iim+i 373 u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij, :,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )374 u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij, :,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )375 u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij, :,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )373 u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) ) 374 u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) ) 375 u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) ) 376 376 ENDDO 377 377 ENDDO -
codes/icosagcm/trunk/src/time/timeloop_gcm.f90
r581 r599 168 168 USE caldyn_mod 169 169 USE advect_tracer_mod 170 USE diagflux_mod 170 171 USE physics_mod 171 172 USE mpipara … … 181 182 REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:) 182 183 183 INTEGER :: ind184 INTEGER :: i t,i,j,l,n, stage185 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time184 REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out 185 INTEGER :: ind, it,i,j,l,n, stage 186 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time 186 187 LOGICAL, PARAMETER :: check_rhodz=.FALSE. 187 188 INTEGER :: start_clock, stop_clock, rate_clock … … 210 211 211 212 IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q) 213 IF(diagflux_on) THEN 214 adv_over_out = itau_adv*(1./itau_out) 215 ELSE 216 adv_over_out = 0. 217 END IF 212 218 213 219 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0) 214 220 215 C ALLtrace_on221 Call trace_on 216 222 217 223 IF (xios_output) THEN ! we must call update_calendar before any XIOS output … … 267 273 CALL HEVI_scheme(it, fluxt_zero) 268 274 END SELECT 269 275 270 276 IF (MOD(it,itau_dissip)==0) THEN 271 277 … … 298 304 299 305 IF(MOD(it,itau_adv)==0) THEN 300 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step 301 fluxt_zero=.TRUE. 302 ! FIXME : check that rhodz is consistent with ps 303 IF (check_rhodz) THEN 306 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step 307 adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt) ! accumulate mass and fluxes if diagflux_on 308 fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme 309 ! At this point advect_tracer has obtained the halos of u and rhodz, 310 ! needed for correct computation of kinetic energy 311 IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt) 312 313 IF (check_rhodz) THEN ! check that rhodz is consistent with ps 304 314 DO ind=1,ndomain 305 315 IF (.NOT. assigned_domain(ind)) CYCLE -
codes/icosagcm/trunk/src/transport/advect.F90
r581 r599 425 425 ! Horizontal transport (S. Dubey, T. Dubos) 426 426 ! Slope-limited van Leer approach with hexagons 427 SUBROUTINE compute_advect_horiz(update_mass, hfluxt,cc,gradq3d, mass,qi)427 SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt) 428 428 USE trace 429 429 USE omp_para 430 430 IMPLICIT NONE 431 LOGICAL, INTENT(IN) :: update_mass 431 LOGICAL, INTENT(IN) :: update_mass, diagflux_on 432 432 REAL(rstd), INTENT(IN) :: gradq3d(iim*jjm,llm,3) 433 433 REAL(rstd), INTENT(IN) :: hfluxt(3*iim*jjm,llm) ! mass flux … … 435 435 REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 436 436 REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm) 437 REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,llm) ! time-integrated tracer flux 437 438 438 439 REAL(rstd) :: dq,dmass,qe,ed(3), newmass … … 441 442 442 443 CALL trace_start("compute_advect_horiz") 443 444 ! evaluate tracer field at point cc using piecewise linear reconstruction 445 ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon 446 ! ne*hfluxt>0 iff outgoing 447 DO l = ll_begin,ll_end 448 449 !DIR$ SIMD 450 DO ij=ij_begin_ext,ij_end_ext 451 452 IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 453 ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 454 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 455 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 456 ELSE 457 ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:) 458 ! qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 459 qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 460 ENDIF 461 qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe 462 463 IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 464 ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 465 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 466 qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 467 ELSE 468 ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:) 469 ! qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 470 qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 471 ENDIF 472 qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 473 474 IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 475 ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 476 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 477 qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 478 ELSE 479 ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:) 480 ! qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 481 qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 482 ENDIF 483 qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe 484 END DO 485 END DO 486 487 ! update q and, if update_mass, update mass 488 DO l = ll_begin,ll_end 489 !DIR$ SIMD 490 DO ij=ij_begin,ij_end 491 ! sign convention as Ringler et al. (2010) eq. 21 492 dmass = hfluxt(ij+u_right,l)*ne(ij,right) & 493 + hfluxt(ij+u_lup,l) *ne(ij,lup) & 494 + hfluxt(ij+u_ldown,l)*ne(ij,ldown) & 495 + hfluxt(ij+u_rup,l) *ne(ij,rup) & 496 + hfluxt(ij+u_left,l) *ne(ij,left) & 497 + hfluxt(ij+u_rdown,l)*ne(ij,rdown) 498 499 dq = qflux(ij+u_right,l) *ne(ij,right) & 500 + qflux(ij+u_lup,l) *ne(ij,lup) & 501 + qflux(ij+u_ldown,l) *ne(ij,ldown) & 502 + qflux(ij+u_rup,l) *ne(ij,rup) & 503 + qflux(ij+u_left,l) *ne(ij,left) & 504 + qflux(ij+u_rdown,l) *ne(ij,rdown) 505 506 507 newmass = mass(ij,l) - dmass/Ai(ij) 508 qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass 509 IF(update_mass) mass(ij,l) = newmass 510 511 END DO 512 END DO 513 444 #include "../kernels/advect_horiz.k90" 514 445 CALL trace_end("compute_advect_horiz") 515 CONTAINS516 517 !====================================================================================518 PURE REAL(rstd) FUNCTION sum2(a1,a2)519 IMPLICIT NONE520 REAL(rstd),INTENT(IN):: a1(3), a2(3)521 sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)522 ! sum2 = 0. ! Godunov scheme523 END FUNCTION sum2524 525 446 END SUBROUTINE compute_advect_horiz 526 447 -
codes/icosagcm/trunk/src/transport/advect_tracer.f90
r548 r599 1 1 MODULE advect_tracer_mod 2 2 USE icosa 3 USE advect_mod 3 4 IMPLICIT NONE 4 5 PRIVATE … … 25 26 26 27 SUBROUTINE init_advect_tracer 27 USE advect_mod28 28 USE omp_para 29 29 REAL(rstd),POINTER :: tangent(:,:) … … 54 54 END SUBROUTINE init_advect_tracer 55 55 56 SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz )57 USE advect_mod58 USE mpipara56 SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz,& 57 frac, f_masst,f_qmasst,f_massfluxt,f_qfluxt) 58 USE omp_para 59 59 USE trace 60 60 USE write_field_mod 61 61 USE tracer_mod 62 IMPLICIT NONE63 64 62 TYPE(t_field),POINTER :: f_hfluxt(:) ! time-integrated horizontal mass flux 65 63 TYPE(t_field),POINTER :: f_wfluxt(:) ! time-integrated vertical mass flux … … 67 65 TYPE(t_field),POINTER :: f_q(:) ! tracer 68 66 TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step 67 REAL(rstd), INTENT(in):: frac ! ratio itau_adv/itau_out or 0. if not diagflux_on 68 TYPE(t_field),POINTER :: f_masst(:) ! time-integrated mass 69 TYPE(t_field),POINTER :: f_qmasst(:) ! time-integrated tracer mass 70 TYPE(t_field),POINTER :: f_massfluxt(:)! time-integrated horizontal mass flux 71 TYPE(t_field),POINTER :: f_qfluxt(:) ! time-integrated horizontal tracer flux 69 72 70 73 REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:) 71 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 74 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:), masst(:,:), qmasst(:,:,:), massfluxt(:,:), qfluxt(:,:,:) 72 75 REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 73 76 ! temporary shared variable for vlz … … 77 80 REAL(rstd),POINTER :: wq(:,:) ! time-integrated flux of q 78 81 79 82 INTEGER :: ind,k, nq_last 80 83 LOGICAL,SAVE :: first=.TRUE. 81 84 !$OMP THREADPRIVATE(first) … … 141 144 CALL send_message(f_cc,req_cc) 142 145 143 144 146 ! horizontal transport - split in two to place transfer of gradq3d 145 147 DO k = 1, nqtot … … 154 156 sqrt_leng=f_sqrt_leng(ind) 155 157 CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v) 156 157 158 END DO 158 159 … … 160 161 CALL wait_message(req_cc) 161 162 CALL wait_message(req_gradq3d) 162 163 163 164 164 DO ind=1,ndomain … … 170 170 rhodz = f_rhodz(ind) 171 171 hfluxt = f_hfluxt(ind) 172 qfluxt = f_qfluxt(ind) 172 173 gradq3d = f_gradq3d(ind) 173 CALL compute_advect_horiz(k==nq_last,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 174 175 IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass 176 qmasst = f_qmasst(ind) 177 qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + & 178 frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k) 179 IF(k==nq_last) THEN 180 masst = f_masst(ind) 181 massfluxt = f_massfluxt(ind) 182 masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end) 183 massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end) 184 END IF 185 END IF 186 CALL compute_advect_horiz(k==nq_last,frac>0., hfluxt,cc,gradq3d, rhodz, q(:,:,k), qfluxt(:,:,k)) 174 187 END DO 188 175 189 ENDIF 176 190 END DO … … 214 228 USE trace 215 229 USE omp_para 216 IMPLICIT NONE217 230 LOGICAL, INTENT(IN) :: update_mass 218 231 REAL(rstd), INTENT(IN) :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux
Note: See TracChangeset
for help on using the changeset viewer.