MODULE caldyn_sw_mod USE icosa USE transfert_mpi_mod, ONLY : t_request PRIVATE TYPE(t_field),POINTER,SAVE :: f_Fe(:) REAL(rstd),POINTER,SAVE :: Fe(:) TYPE(t_field),POINTER,SAVE :: f_Ki(:) REAL(rstd),POINTER,SAVE :: Ki(:) TYPE(t_field),POINTER,SAVE :: f_qv(:) REAL(rstd),POINTER,SAVE :: qv(:) TYPE(t_field),POINTER,SAVE :: f_qv2(:) REAL(rstd),POINTER,SAVE :: qv2(:) TYPE(t_field),POINTER,SAVE :: f_t_tmp(:) REAL(rstd),POINTER,SAVE :: t_tmp(:) TYPE(t_field),POINTER,SAVE :: f_u_tmp(:) REAL(rstd),POINTER,SAVE :: u_tmp(:) TYPE(t_field),POINTER,SAVE :: f_z_tmp(:) REAL(rstd),POINTER,SAVE :: z_tmp(:) TYPE(t_field),POINTER,SAVE :: f_dZ(:) REAL(rstd),POINTER,SAVE :: dZ(:) TYPE(t_field),POINTER,SAVE :: f_diff_dhv(:) REAL(rstd),POINTER,SAVE :: diff_dhv(:) TYPE(t_request),POINTER :: req(:) TYPE(t_request),POINTER :: req_u(:) PUBLIC :: allocate_caldyn,swap_caldyn,caldyn,write_caldyn,Compute_PV,Compute_enstrophy CONTAINS SUBROUTINE allocate_caldyn USE icosa IMPLICIT NONE INTEGER :: ind,i,j CALL allocate_field(f_Fe,field_u,type_real) CALL allocate_field(f_Ki,field_t,type_real) CALL allocate_field(f_qv,field_z,type_real) CALL allocate_field(f_qv2,field_z,type_real) CALL allocate_field(f_t_tmp,field_t,type_real) CALL allocate_field(f_u_tmp,field_u,type_real) CALL allocate_field(f_z_tmp,field_z,type_real) CALL allocate_field(f_dZ,field_z,type_real) CALL allocate_field(f_diff_dhv,field_z,type_real) CALL create_request(field_t,req) DO ind=1,ndomain DO i=ii_begin,ii_end+1 CALL request_add_point(ind,i,jj_begin-1,req) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_end+1,j,req) ENDDO DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_end+1,req) ENDDO DO j=jj_begin,jj_end+1 CALL request_add_point(ind,ii_begin-1,j,req) ENDDO DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_begin,req) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_end,j,req) ENDDO DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_end,req) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_begin,j,req) ENDDO ENDDO CALL create_request(field_u,req_u) DO ind=1,ndomain DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_begin-1,req_u,rup) CALL request_add_point(ind,i+1,jj_begin-1,req_u,lup) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_end+1,j,req_u,left) CALL request_add_point(ind,ii_end+1,j-1,req_u,lup) ENDDO DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_end+1,req_u,ldown) CALL request_add_point(ind,i-1,jj_end+1,req_u,rdown) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_begin-1,j,req_u,right) CALL request_add_point(ind,ii_begin-1,j+1,req_u,rdown) ENDDO ENDDO END SUBROUTINE allocate_caldyn SUBROUTINE swap_caldyn(ind) USE icosa IMPLICIT NONE INTEGER,INTENT(IN) :: ind Fe=f_Fe(ind) Ki=f_Ki(ind) qv=f_qv(ind) qv2=f_qv2(ind) t_tmp=f_t_tmp(ind) u_tmp=f_u_tmp(ind) z_tmp=f_z_tmp(ind) dZ=f_dZ(ind) diff_dhv=f_diff_dhv(ind) END SUBROUTINE swap_caldyn SUBROUTINE caldyn(f_h, f_u, f_dh, f_du) USE icosa IMPLICIT NONE TYPE(t_field),POINTER :: f_h(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_dh(:) TYPE(t_field),POINTER :: f_du(:) REAL(rstd),POINTER :: h(:) REAL(rstd),POINTER :: u(:) REAL(rstd),POINTER :: dh(:) REAL(rstd),POINTER :: du(:) INTEGER :: ind INTEGER,SAVE :: it=0 CALL transfert_request(f_h,req_i1) CALL transfert_request(f_u,req_e1) CALL transfert_request(f_u,req_e1) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) CALL swap_caldyn(ind) h=f_h(ind) u=f_u(ind) dh=f_dh(ind) du=f_du(ind) CALL compute_caldyn(h, u, dh, du) ENDDO IF (mod(it,240)==0) THEN CALL writefield("h",f_h) CALL writefield("dh",f_dh) CALL Compute_enstrophy ! CALL Compute_PV ENDIF it=it+1 END SUBROUTINE caldyn SUBROUTINE compute_caldyn(hi,ue,dhi,due) USE icosa IMPLICIT NONE REAL(rstd),INTENT(IN) :: hi(iim*jjm) REAL(rstd),INTENT(IN) :: ue(iim*3*jjm) REAL(rstd),INTENT(OUT):: dhi(iim*jjm) REAL(rstd),INTENT(OUT):: due(iim*3*jjm) REAL(rstd) :: Fep(iim*3*jjm) INTEGER :: i,j,n REAL(rstd) :: etav,hv DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 n=(j-1)*iim+i Fe(n+u_right)=0.5*(hi(n)+hi(n+t_right))*ue(n+u_right) Fe(n+u_lup)=0.5*(hi(n)+hi(n+t_lup))*ue(n+u_lup) Fe(n+u_ldown)=0.5*(hi(n)+hi(n+t_ldown))*ue(n+u_ldown) u_tmp(n+u_right)=ue(n+u_right) u_tmp(n+u_lup)=ue(n+u_lup) u_tmp(n+u_ldown)=ue(n+u_ldown) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i dhi(n)=-1./Ai(n)*(ne(n,right)*Fe(n+u_right)*le(n+u_right) + & ne(n,rup)*Fe(n+u_rup)*le(n+u_rup) + & ne(n,lup)*Fe(n+u_lup)*le(n+u_lup) + & ne(n,left)*Fe(n+u_left)*le(n+u_left) + & ne(n,ldown)*Fe(n+u_ldown)*le(n+u_ldown) + & ne(n,rdown)*Fe(n+u_rdown)*le(n+u_rdown)) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i Ki(n)=1/(4*Ai(n))*(le(n+u_right)*de(n+u_right)*ue(n+u_right)**2 + & le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)**2 + & le(n+u_lup)*de(n+u_lup)*ue(n+u_lup)**2 + & le(n+u_left)*de(n+u_left)*ue(n+u_left)**2 + & le(n+u_ldown)*de(n+u_ldown)*ue(n+u_ldown)**2 + & le(n+u_rdown)*de(n+u_rdown)*ue(n+u_rdown)**2 ) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i due(n+u_right)=1/de(n+u_right)*(ne(n,right)*(g*(hi(n)+bi(n))+Ki(n))+ & ne(n+t_right,left)*(g*(hi(n+t_right)+bi(n+t_right))+Ki(n+t_right)) ) due(n+u_lup)=1/de(n+u_lup)*(ne(n,lup)*(g*(hi(n)+bi(n))+Ki(n))+ & ne(n+t_lup,rdown)*(g*(hi(n+t_lup)+bi(n+t_lup))+Ki(n+t_lup)) ) due(n+u_ldown)=1/de(n+u_ldown)*(ne(n,ldown)*(g*(hi(n)+bi(n))+Ki(n))+ & ne(n+t_ldown,rup)*(g*(hi(n+t_ldown)+bi(n+t_ldown))+Ki(n+t_ldown)) ) ENDDO ENDDO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 n=(j-1)*iim+i etav= 1./Av(n+z_up)*( ne(n,rup)*ue(n+u_rup)*de(n+u_rup) & + ne(n+t_rup,left)*ue(n+t_rup+u_left)*de(n+t_rup+u_left) & - ne(n,lup)*ue(n+u_lup)*de(n+u_lup) ) hv= Riv2(n,vup)*hi(n)+ & Riv2(n+t_rup,vldown)*hi(n+t_rup)+ & Riv2(n+t_lup,vrdown)*hi(n+t_lup) qv(n+z_up)=(etav+fv(n+z_up))/hv qv2(n+z_up)= (etav+fv(n+z_down))**2/hv z_tmp(n+z_up)= Av(n+z_up) etav = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown)*de(n+u_ldown) & + ne(n+t_ldown,right)*ue(n+t_ldown+u_right)*de(n+t_ldown+u_right) & - ne(n,rdown)*ue(n+u_rdown)*de(n+u_rdown) ) hv = Riv2(n,vdown)*hi(n)+ & Riv2(n+t_ldown,vrup)*hi(n+t_ldown)+ & Riv2(n+t_rdown,vlup)*hi(n+t_rdown) qv(n+z_down)=(etav+fv(n+z_down))/hv qv2(n+z_down)= (etav+fv(n+z_down))**2/hv z_tmp(n+z_down)=Av(n+z_down) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i due(n+u_right) = due(n+u_right)+0.5*(qv(n+z_rdown)+qv(n+z_rup))/de(n+u_right) * & ( wee(n+u_right,1,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_right,2,1)*le(n+u_lup)*Fe(n+u_lup)+ & wee(n+u_right,3,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_right,4,1)*le(n+u_ldown)*Fe(n+u_ldown)+ & wee(n+u_right,5,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_right,1,2)*le(n+t_right+u_ldown)*Fe(n+t_right+u_ldown)+ & wee(n+u_right,2,2)*le(n+t_right+u_rdown)*Fe(n+t_right+u_rdown)+ & wee(n+u_right,3,2)*le(n+t_right+u_right)*Fe(n+t_right+u_right)+ & wee(n+u_right,4,2)*le(n+t_right+u_rup)*Fe(n+t_right+u_rup)+ & wee(n+u_right,5,2)*le(n+t_right+u_lup)*Fe(n+t_right+u_lup) ) due(n+u_lup) = due(n+u_lup) + 0.5*(qv(n+z_up)+qv(n+z_lup))/de(n+u_lup) * & ( wee(n+u_lup,1,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_lup,2,1)*le(n+u_ldown)*Fe(n+u_ldown)+ & wee(n+u_lup,3,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_lup,4,1)*le(n+u_right)*Fe(n+u_right)+ & wee(n+u_lup,5,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_lup,1,2)*le(n+t_lup+u_right)*Fe(n+t_lup+u_right)+ & wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*Fe(n+t_lup+u_rup)+ & wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*Fe(n+t_lup+u_lup)+ & wee(n+u_lup,4,2)*le(n+t_lup+u_left)*Fe(n+t_lup+u_left)+ & wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*Fe(n+t_lup+u_ldown) ) due(n+u_ldown) = due(n+u_ldown)+ 0.5*(qv(n+z_ldown)+qv(n+z_down))/de(n+u_ldown) * & ( wee(n+u_ldown,1,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_ldown,2,1)*le(n+u_right)*Fe(n+u_right)+ & wee(n+u_ldown,3,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_ldown,4,1)*le(n+u_lup)*Fe(n+u_lup)+ & wee(n+u_ldown,5,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*Fe(n+t_ldown+u_lup)+ & wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*Fe(n+t_ldown+u_left)+ & wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*Fe(n+t_ldown+u_ldown)+ & wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*Fe(n+t_ldown+u_rdown)+ & wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*Fe(n+t_ldown+u_right) ) ENDDO ENDDO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 n=(j-1)*iim+i dZ(n+z_up) = 2*qv(n+z_up)/Av(n+z_up)*( ne(n,rup)*due(n+u_rup)*de(n+u_rup) & + ne(n+t_rup,left)*due(n+t_rup+u_left)*de(n+t_rup+u_left) & - ne(n,lup)*due(n+u_lup)*de(n+u_lup) ) dZ(n+z_up)= dZ(n+z_up)-(Riv2(n,vup)*dhi(n)+ & Riv2(n+t_rup,vldown)*dhi(n+t_rup)+ & Riv2(n+t_lup,vrdown)*dhi(n+t_lup))*qv(n+z_up)**2 dZ(n+z_down) = 2*qv(n+z_down)/Av(n+z_down)*( ne(n,ldown)*due(n+u_ldown)*de(n+u_ldown) & + ne(n+t_ldown,right)*due(n+t_ldown+u_right)*de(n+t_ldown+u_right) & - ne(n,rdown)*due(n+u_rdown)*de(n+u_rdown) ) dZ(n+z_down)= dZ(n+z_down)- (Riv2(n,vdown)*dhi(n)+ & Riv2(n+t_ldown,vrup)*dhi(n+t_ldown)+ & Riv2(n+t_rdown,vlup)*dhi(n+t_rdown))*qv(n+z_down)**2 ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i Fep(n+u_right) = 1./de(n+u_right) * & ( wee(n+u_right,1,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_right,2,1)*le(n+u_lup)*Fe(n+u_lup)+ & wee(n+u_right,3,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_right,4,1)*le(n+u_ldown)*Fe(n+u_ldown)+ & wee(n+u_right,5,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_right,1,2)*le(n+t_right+u_ldown)*Fe(n+t_right+u_ldown)+ & wee(n+u_right,2,2)*le(n+t_right+u_rdown)*Fe(n+t_right+u_rdown)+ & wee(n+u_right,3,2)*le(n+t_right+u_right)*Fe(n+t_right+u_right)+ & wee(n+u_right,4,2)*le(n+t_right+u_rup)*Fe(n+t_right+u_rup)+ & wee(n+u_right,5,2)*le(n+t_right+u_lup)*Fe(n+t_right+u_lup) ) Fep(n+u_lup) = 1/de(n+u_lup) * & ( wee(n+u_lup,1,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_lup,2,1)*le(n+u_ldown)*Fe(n+u_ldown)+ & wee(n+u_lup,3,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_lup,4,1)*le(n+u_right)*Fe(n+u_right)+ & wee(n+u_lup,5,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_lup,1,2)*le(n+t_lup+u_right)*Fe(n+t_lup+u_right)+ & wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*Fe(n+t_lup+u_rup)+ & wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*Fe(n+t_lup+u_lup)+ & wee(n+u_lup,4,2)*le(n+t_lup+u_left)*Fe(n+t_lup+u_left)+ & wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*Fe(n+t_lup+u_ldown) ) Fep(n+u_ldown) = 1./de(n+u_ldown) * & ( wee(n+u_ldown,1,1)*le(n+u_rdown)*Fe(n+u_rdown)+ & wee(n+u_ldown,2,1)*le(n+u_right)*Fe(n+u_right)+ & wee(n+u_ldown,3,1)*le(n+u_rup)*Fe(n+u_rup)+ & wee(n+u_ldown,4,1)*le(n+u_lup)*Fe(n+u_lup)+ & wee(n+u_ldown,5,1)*le(n+u_left)*Fe(n+u_left)+ & wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*Fe(n+t_ldown+u_lup)+ & wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*Fe(n+t_ldown+u_left)+ & wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*Fe(n+t_ldown+u_ldown)+ & wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*Fe(n+t_ldown+u_rdown)+ & wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*Fe(n+t_ldown+u_right) ) ENDDO ENDDO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 n=(j-1)*iim+i diff_dhv(n+z_up)= 1./Av(n+z_up)*( ne(n,rup)*Fep(n+u_rup)*de(n+u_rup) & + ne(n+t_rup,left)*Fep(n+t_rup+u_left)*de(n+t_rup+u_left) & - ne(n,lup)*Fep(n+u_lup)*de(n+u_lup) ) diff_dhv(n+z_up)=diff_dhv(n+z_up)-(Riv2(n,vup)*dhi(n)+ & Riv2(n+t_rup,vldown)*dhi(n+t_rup)+ & Riv2(n+t_lup,vrdown)*dhi(n+t_lup)) diff_dhv(n+z_down) = 1./Av(n+z_down)*( ne(n,ldown)*Fep(n+u_ldown)*de(n+u_ldown) & + ne(n+t_ldown,right)*Fep(n+t_ldown+u_right)*de(n+t_ldown+u_right) & - ne(n,rdown)*Fep(n+u_rdown)*de(n+u_rdown) ) diff_dhv(n+z_down) = diff_dhv(n+z_down) - (Riv2(n,vdown)*dhi(n)+ & Riv2(n+t_ldown,vrup)*dhi(n+t_ldown)+ & Riv2(n+t_rdown,vlup)*dhi(n+t_rdown)) ENDDO ENDDO END SUBROUTINE compute_caldyn SUBROUTINE write_caldyn USE icosa IMPLICIT NONE INTEGER :: ind,i,j,n CALL transfert_request(f_u_tmp,req_u) DO ind=1,ndomain CALL swap_caldyn(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i Ki(n)=1/(4*Ai(n))*(le(n+u_right)*de(n+u_right)*u_tmp(n+u_right)**2 + & le(n+u_rup)*de(n+u_rup)*u_tmp(n+u_rup)**2 + & le(n+u_lup)*de(n+u_lup)*u_tmp(n+u_lup)**2 + & le(n+u_left)*de(n+u_left)*u_tmp(n+u_left)**2 + & le(n+u_ldown)*de(n+u_ldown)*u_tmp(n+u_ldown)**2 + & le(n+u_rdown)*de(n+u_rdown)*u_tmp(n+u_rdown)**2 ) ENDDO ENDDO ENDDO CALL writefield("Ki",f_Ki) END SUBROUTINE write_caldyn SUBROUTINE Compute_PV USE icosa IMPLICIT NONE REAL(rstd) :: PV INTEGER :: i,j,n,ind PV=0 DO ind=1,ndomain CALL swap_caldyn(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=(j-1)*iim+i qv(n+z_down)=qv(n+z_down)*Av(n+z_down) PV=PV+qv(n+z_down) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=(j-1)*iim+i qv(n+z_up)=qv(n+z_up)*Av(n+z_up) PV=PV+qv(n+z_up) ENDDO ENDDO ENDDO CALL writefield("PV",f_qv) PRINT *,"PV=",PV END SUBROUTINE Compute_PV SUBROUTINE Compute_enstrophy USE icosa IMPLICIT NONE REAL(rstd) :: Es INTEGER :: i,j,n,ind Es=0 DO ind=1,ndomain CALL swap_caldyn(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=(j-1)*iim+i Es=Es+qv2(n+z_down)*Av(n+z_down) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=(j-1)*iim+i Es=Es+qv2(n+z_up)*Av(n+z_up) ENDDO ENDDO ENDDO PRINT *,"Enstrophy=",Es CALL writefield("Ens",f_qv2) Es=0 DO ind=1,ndomain CALL swap_caldyn(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=(j-1)*iim+i Es=Es+dZ(n+z_down)*Av(n+z_down) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=(j-1)*iim+i Es=Es+dZ(n+z_up)*Av(n+z_up) ENDDO ENDDO ENDDO PRINT *,"D/dt Enstrophy=",Es CALL writefield("dEns",f_dZ) CALL writefield("diff_dhv",f_diff_dhv) END SUBROUTINE Compute_enstrophy END MODULE caldyn_sw_mod