MODULE dissip_mod USE genmod USE field_mod USE transfert_mod TYPE(t_field),POINTER,SAVE :: f_gradrot(:) TYPE(t_request),POINTER :: req_dissip(:) INTEGER,PARAMETER :: nitergdiv=1 INTEGER,PARAMETER :: nitergrot=1 REAL :: tetagdiv REAL :: tetagrot REAL :: cdivu REAL :: crot INTEGER :: idissip REAL :: dtdissip CONTAINS SUBROUTINE allocate_dissip IMPLICIT NONE CALL allocate_field(f_gradrot,field_u,type_real) END SUBROUTINE allocate_dissip SUBROUTINE init_dissip(dt) USE domain_mod USE dimensions USE geometry USE metric USE ioipsl IMPLICIT NONE REAL,INTENT(IN) :: dt TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_du(:) REAL(rstd),POINTER :: u(:) REAL(rstd),POINTER :: du(:) REAL(rstd) :: dumax,dumaxmax,dumaxm1 REAL(rstd) :: r INTEGER :: i,j,n,ind,it tetagdiv=5000 CALL getin("tetagdiv",tetagdiv) tetagrot=5000 CALL getin("tetagrot",tetagrot) CALL allocate_dissip CALL allocate_field(f_u,field_u,type_real) CALL allocate_field(f_du,field_u,type_real) CALL create_request(field_u,req_dissip) DO ind=1,ndomain DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup) CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_end+1,j,req_dissip,left) CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup) ENDDO DO i=ii_begin,ii_end CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown) CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown) ENDDO DO j=jj_begin,jj_end CALL request_add_point(ind,ii_begin-1,j,req_dissip,right) CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown) ENDDO DO i=ii_begin+1,ii_end-1 CALL request_add_point(ind,i,jj_begin,req_dissip,right) CALL request_add_point(ind,i,jj_end,req_dissip,right) ENDDO DO j=jj_begin+1,jj_end-1 CALL request_add_point(ind,ii_begin,j,req_dissip,rup) CALL request_add_point(ind,ii_end,j,req_dissip,rup) ENDDO CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left) CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown) CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left) CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown) ENDDO cdivu=-1 crot=-1 CALL RANDOM_SEED() DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) u=f_u(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i CALL RANDOM_NUMBER(r) u(n+u_right)=r-0.5 CALL RANDOM_NUMBER(r) u(n+u_lup)=r-0.5 CALL RANDOM_NUMBER(r) u(n+u_ldown)=r-0.5 ENDDO ENDDO ENDDO CALL transfert_request(f_u,req_dissip) CALL transfert_request(f_u,req_dissip) DO it=1,500 dumaxm1=dumax dumax=0 DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) u=f_u(ind) du=f_du(ind) CALL gradiv(u,du) ENDDO CALL transfert_request(f_du,req_dissip) CALL transfert_request(f_du,req_dissip) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) du=f_du(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right))) if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup))) if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown))) ENDDO ENDDO ENDDO u=du/dumax PRINT *,"gradiv : it :",it ,": dumax",(dumax+dumaxm1)/2 ENDDO PRINT *,"gradiv : dumax",(dumax+dumaxm1)/2 cdivu=dumax**(-1./nitergdiv) PRINT *,"cdivu : ",cdivu DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) u=f_u(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i CALL RANDOM_NUMBER(r) u(n+u_right)=r-0.5 CALL RANDOM_NUMBER(r) u(n+u_lup)=r-0.5 CALL RANDOM_NUMBER(r) u(n+u_ldown)=r-0.5 ENDDO ENDDO ENDDO CALL transfert_request(f_u,req_dissip) CALL transfert_request(f_u,req_dissip) DO it=1,500 dumaxm1=dumax dumax=0 DO ind=1,ndomain u=f_u(ind) du=f_du(ind) CALL gradrot(u,du) ENDDO CALL transfert_request(f_du,req_dissip) CALL transfert_request(f_du,req_dissip) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) du=f_du(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right))) if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup))) if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown))) ENDDO ENDDO ENDDO u=du/dumax PRINT *,"gradrot : it :",it ,": dumax",(dumax+dumaxm1)/2 ENDDO PRINT *,"gradrot : dumax",(dumax+dumaxm1)/2 crot=dumax**(-1/nitergrot) PRINT *,"crot : ",crot idissip=INT(MIN(tetagdiv,tetagrot)/(2.*dt)) idissip=MAX(1,idissip) dtdissip=idissip*dt PRINT *,"idissip",idissip," dtdissip ",dtdissip END SUBROUTINE init_dissip SUBROUTINE dissip(f_ue,f_due) USE domain_mod USE dimensions USE geometry USE metric IMPLICIT NONE TYPE(t_field),POINTER :: f_ue(:) TYPE(t_field),POINTER :: f_due(:) REAL(rstd),POINTER :: ue(:) REAL(rstd),POINTER :: due(:) REAL(rstd),POINTER :: gradrot_e(:) INTEGER :: ind REAL :: v CALL transfert_request(f_ue,req_dissip) CALL transfert_request(f_ue,req_dissip) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ue=f_ue(ind) due=f_due(ind) gradrot_e=f_gradrot(ind) CALL gradiv(ue,due) due=due*dtdissip/tetagdiv CALL gradrot(ue,gradrot_e) due=due+gradrot_e*dtdissip/tetagrot ENDDO END SUBROUTINE dissip SUBROUTINE gradiv(ue,gradivu_e) USE domain_mod USE dimensions USE geometry USE metric IMPLICIT NONE REAL(rstd),INTENT(IN) :: ue(iim*3*jjm) REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm) REAL(rstd) :: divu_i(iim*jjm) INTEGER :: i,j,n DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i divu_i(n)=-1./Ai(n)*(ne(n,right)*ue(n+u_right)*le(n+u_right) + & ne(n,rup)*ue(n+u_rup)*le(n+u_rup) + & ne(n,lup)*ue(n+u_lup)*le(n+u_lup) + & ne(n,left)*ue(n+u_left)*le(n+u_left) + & ne(n,ldown)*ue(n+u_ldown)*le(n+u_ldown) + & ne(n,rdown)*ue(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 gradivu_e(n+u_right)=1/de(n+u_right)*(ne(n,right)*divu_i(n)+ ne(n+t_right,left)*divu_i(n+t_right) ) gradivu_e(n+u_lup)=1/de(n+u_lup)*(ne(n,lup)*divu_i(n)+ ne(n+t_lup,rdown)*divu_i(n+t_lup )) gradivu_e(n+u_ldown)=1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n)+ne(n+t_ldown,rup)*divu_i(n+t_ldown) ) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i gradivu_e(n+u_right)=-gradivu_e(n+u_right)*cdivu gradivu_e(n+u_lup)=-gradivu_e(n+u_lup)*cdivu gradivu_e(n+u_ldown)=-gradivu_e(n+u_ldown)*cdivu ENDDO ENDDO END SUBROUTINE gradiv SUBROUTINE gradrot(ue,gradrot_e) USE domain_mod USE dimensions USE geometry USE metric IMPLICIT NONE REAL(rstd),INTENT(IN) :: ue(iim*3*jjm) REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm) REAL(rstd) :: rot_v(2*iim*jjm) INTEGER :: i,j,n DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 n=(j-1)*iim+i rot_v(n+z_up)= 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) ) rot_v(n+z_down) = 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) ) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i gradrot_e(n+u_right)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown)-rot_v(n+z_rup)) gradrot_e(n+u_lup)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up)-rot_v(n+z_lup)) gradrot_e(n+u_ldown)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown)-rot_v(n+z_down)) ENDDO ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end n=(j-1)*iim+i gradrot_e(n+u_right)=-gradrot_e(n+u_right)*crot gradrot_e(n+u_lup)=-gradrot_e(n+u_lup)*crot gradrot_e(n+u_ldown)=-gradrot_e(n+u_ldown)*crot ENDDO ENDDO END SUBROUTINE END MODULE dissip_mod