MODULE caldyn_adv_mod USE icosa TYPE(t_field),POINTER :: f_out(:) REAL(rstd),POINTER :: out(:,:) TYPE(t_field),POINTER :: f_out_u(:) REAL(rstd),POINTER :: out_u(:,:) TYPE(t_field),POINTER :: f_out_z(:) REAL(rstd),POINTER :: out_z(:,:) INTEGER :: itau_out CONTAINS SUBROUTINE init_caldyn(dt) USE icosa IMPLICIT NONE REAL(rstd),INTENT(IN) :: dt INTEGER :: write_period CALL allocate_caldyn CALL getin('write_period',write_period) itau_out=INT(write_period/dt) CALL allocate_caldyn END SUBROUTINE init_caldyn SUBROUTINE allocate_caldyn USE icosa IMPLICIT NONE CALL allocate_field(f_out,field_t,type_real,llm) CALL allocate_field(f_out_u,field_u,type_real,llm) CALL allocate_field(f_out_z,field_z,type_real,llm) END SUBROUTINE allocate_caldyn SUBROUTINE swap_caldyn(ind) IMPLICIT NONE INTEGER,INTENT(IN) :: ind out=f_out(ind) out_u=f_out_u(ind) out_z=f_out_z(ind) END SUBROUTINE swap_caldyn SUBROUTINE check_mass_conservation(f_ps,f_dps) USE icosa IMPLICIT NONE TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_dps(:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: dps(:) REAL(rstd) :: mass_tot,dmass_tot INTEGER :: ind,i,j,ij mass_tot=0 dmass_tot=0 CALL transfert_request(f_dps,req_i1) CALL transfert_request(f_ps,req_i1) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) dps=f_dps(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i IF (domain(ind)%own(i,j)) THEN mass_tot=mass_tot+ps(ij)*Ai(ij)/g dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g ENDIF ENDDO ENDDO ENDDO PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot END SUBROUTINE check_mass_conservation SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) USE icosa USE vorticity_mod USE kinetic_mod USE theta2theta_rhodz_mod IMPLICIT NONE INTEGER,INTENT(IN) :: it TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_dps(:) TYPE(t_field),POINTER :: f_dtheta_rhodz(:) TYPE(t_field),POINTER :: f_du(:) REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: theta_rhodz(:,:) REAL(rstd),POINTER :: u(:,:) REAL(rstd),POINTER :: dps(:) REAL(rstd),POINTER :: dtheta_rhodz(:,:) REAL(rstd),POINTER :: du(:,:) INTEGER :: ind CALL transfert_request(f_phis,req_i1) CALL transfert_request(f_ps,req_i1) CALL transfert_request(f_theta_rhodz,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) phis=f_phis(ind) ps=f_ps(ind) theta_rhodz=f_theta_rhodz(ind) u=f_u(ind) dps=f_dps(ind) dtheta_rhodz=f_dtheta_rhodz(ind) du=f_du(ind) !$OMP PARALLEL DEFAULT(SHARED) CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du) !$OMP END PARALLEL ENDDO CALL transfert_request(f_out_u,req_e1) CALL transfert_request(f_out_u,req_e1) IF (mod(it,itau_out)==0 ) THEN CALL writefield("ps",f_ps) ENDIF ! CALL check_mass_conservation(f_ps,f_dps) END SUBROUTINE caldyn SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du) USE icosa USE disvert_mod IMPLICIT NONE REAL(rstd),INTENT(IN) :: phis(iim*jjm) REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) REAL(rstd),INTENT(IN) :: ps(iim*jjm) REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm) REAL(rstd),INTENT(OUT):: dps(iim*jjm) INTEGER :: i,j,ij,l REAL(rstd) :: ww,uu REAL(rstd) :: delta REAL(rstd) :: etav,hv REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:) ! potential temperature REAL(rstd),ALLOCATABLE,SAVE :: p(:,:) ! pression REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:) ! Exner function REAL(rstd),ALLOCATABLE,SAVE :: pks(:) ! Intermediate variable to compute exner function REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:) REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:) ! REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:) ! geopotential REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:) ! mass REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:) ! mass density REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:) ! mass flux REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:) ! mass flux convergence REAL(rstd),ALLOCATABLE,SAVE :: w(:,:) ! vertical velocity REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:) ! potential velocity REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:) ! bernouilli term LOGICAL,SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) !$OMP BARRIER !$OMP MASTER IF (first) THEN ALLOCATE(theta(iim*jjm,llm)) ! potential temperature ALLOCATE(p(iim*jjm,llm+1)) ! pression ALLOCATE(pk(iim*jjm,llm)) ! Exner function ALLOCATE(pks(iim*jjm)) ALLOCATE(alpha(iim*jjm,llm)) ALLOCATE(beta(iim*jjm,llm)) ALLOCATE(phi(iim*jjm,llm)) ! geopotential ALLOCATE(mass(iim*jjm,llm)) ! mass ALLOCATE(rhodz(iim*jjm,llm)) ! mass density ALLOCATE(Fe(3*iim*jjm,llm)) ! mass flux ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux ALLOCATE(convm(iim*jjm,llm)) ! mass flux convergence ALLOCATE(w(iim*jjm,llm)) ! vertical velocity ALLOCATE(qv(2*iim*jjm,llm)) ! potential velocity ALLOCATE(berni(iim*jjm,llm)) ! bernouilli term first=.FALSE. ENDIF !$OMP END MASTER !$OMP BARRIER !!! Compute pression dtheta_rhodz(:,:)=0. du(:,:)=0. DO l = 1, llm+1 !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i p(ij,l) = ap(l) + bp(l) * ps(ij) ENDDO ENDDO ENDDO !!! Compute mass DO l = 1, llm !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g rhodz(ij,l) = mass(ij,l) / Ai(ij) ENDDO ENDDO ENDDO !!! Compute mass flux !! question à thomas : meilleure pondération de la masse sur les liens ? DO l = 1, llm !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) ENDDO ENDDO ENDDO !!! mass flux convergence computation ! horizontal convergence DO l = 1, llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i !signe ? convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l) + & ne(ij,rup)*Fe(ij+u_rup,l) + & ne(ij,lup)*Fe(ij+u_lup,l) + & ne(ij,left)*Fe(ij+u_left,l) + & ne(ij,ldown)*Fe(ij+u_ldown,l) + & ne(ij,rdown)*Fe(ij+u_rdown,l)) ENDDO ENDDO ENDDO ! vertical integration from up to down DO l = llm-1, 1, -1 !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i convm(ij,l) = convm(ij,l) + convm(ij,l+1) ENDDO ENDDO ENDDO !!! Compute dps !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i dps(ij)=-convm(ij,1) * g ENDDO ENDDO END SUBROUTINE compute_caldyn END MODULE caldyn_adv_mod