MODULE caldyn_gcm_mod USE genmod USE field_mod USE transfert_mod 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(:,:) CONTAINS SUBROUTINE allocate_caldyn USE domain_mod USE dimensions USE geometry USE metric IMPLICIT NONE INTEGER :: ind,i,j 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 domain_mod USE dimensions USE geometry 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(f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) USE domain_mod USE dimensions USE grid_param USE geometry USE metric USE write_field USE vorticity_mod USE kinetic_mod IMPLICIT NONE 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 INTEGER,SAVE :: it=0 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) ! CALL vorticity(f_u,f_out_z) 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) ! CALL vorticity(f_u,f_out_z) ! CALL kinetic(f_du,f_out) IF (mod(it,72)==0 ) THEN CALL writefield("ps",f_ps) CALL writefield("dps",f_dps) ! CALL writefield("theta_rhodz",f_theta_rhodz) CALL writefield("dtheta_rhodz",f_dtheta_rhodz) ! CALL writefield("vort",f_out_z) CALL writefield("theta",f_out) ! CALL writefield("out",f_out) ! DO ind=1,ndomain ! CALL writefield("Ki",f_out,ind) ! CALL writefield("vort",f_out_z,ind) ! ENDDO ENDIF ! CALL check_mass_conservation(f_ps,f_dps) it=it+1 END SUBROUTINE caldyn SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du) USE domain_mod USE dimensions USE grid_param USE geometry USE metric 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) :: theta(iim*jjm,llm) ! potential temperature ! REAL(rstd) :: p(iim*jjm,llm+1) ! pression ! REAL(rstd) :: pk(iim*jjm,llm) ! Exner function ! REAL(rstd) :: pks(iim*jjm) !! Intermediate variable to compute exner function ! REAL(rstd) :: alpha(iim*jjm,llm) ! REAL(rstd) :: beta(iim*jjm,llm) !! ! REAL(rstd) :: phi(iim*jjm,llm) ! geopotential ! REAL(rstd) :: mass(iim*jjm,llm) ! mass ! REAL(rstd) :: rhodz(iim*jjm,llm) ! mass density ! REAL(rstd) :: Fe(3*iim*jjm,llm) ! mass flux ! REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux ! REAL(rstd) :: convm(iim*jjm,llm) ! mass flux convergence ! REAL(rstd) :: w(iim*jjm,llm) ! vertical velocity ! REAL(rstd) :: qv(2*iim*jjm,llm) ! potential velocity ! REAL(rstd) :: berni(iim*jjm,llm) ! bernouilli term 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 ! du(:,:)=0 ! theta=1e10 ! p=1e10 ! pk=1e10 ! pks=1e10 ! alpha=1e10 ! beta=1e10 ! phi=1e10 ! mass=1e10 ! rhodz=1e10 ! Fe=1e10 ! Ftheta=1e10 ! convm=1e10 ! w=1e10 ! qv=1e10 ! berni=1e10 !!! Compute pression 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 Exnher function !! Compute Alpha and Beta ! for llm layer !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i alpha(ij,llm) = 0. beta (ij,llm) = 1./ (1+ 2*kappa) ENDDO ENDDO ! for other layer DO l = llm-1 , 2 , -1 !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) ) alpha(ij,l) = - p(ij,l+1) / delta * alpha(ij,l+1) beta (ij,l) = p(ij,l ) / delta ENDDO ENDDO ENDDO !! Compute pk ! for first layer !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i pks(ij) = cpp * ( ps(ij)/preff ) ** kappa pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) / & ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) ) ENDDO ENDDO ! for other layers DO l = 2, llm !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 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 theta 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 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) ENDDO ENDDO ENDDO !!! Compute geopotential ! for first layer !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) ENDDO ENDDO ! for other layers DO l = 2, llm !$OMP DO DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l) + theta(ij,l-1) ) & * ( pk(ij,l-1) - pk(ij,l) ) 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 !!! fisrt composante dtheta ! Flux on the edge 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 Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l) Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l) Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l) ENDDO ENDDO ENDDO ! compute divergence DO l = 1, llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ! signe ? attention d (rho theta dz) dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l) + & ne(ij,rup)*Ftheta(ij+u_rup,l) + & ne(ij,lup)*Ftheta(ij+u_lup,l) + & ne(ij,left)*Ftheta(ij+u_left,l) + & ne(ij,ldown)*Ftheta(ij+u_ldown,l) + & ne(ij,rdown)*Ftheta(ij+u_rdown,l)) 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 DO l = 1, llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i out(ij,l)=theta(ij,l)-288 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 !!! Compute vertical velocity DO l = 1,llm-1 !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) ENDDO ENDDO ENDDO !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i w(ij,1) = 0. ENDDO ENDDO !!! Compute potential vorticity 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 etav= 1./Av(ij+z_up)*( ne(ij,rup) * u(ij+u_rup,l) * de(ij+u_rup) & + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) & - ne(ij,lup) * u(ij+u_lup,l) * de(ij+u_lup) ) hv = Riv2(ij,vup) * rhodz(ij,l) & + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv etav = 1./Av(ij+z_down)*( ne(ij,ldown) * u(ij+u_ldown,l) * de(ij+u_ldown) & + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) & - ne(ij,rdown) * u(ij+u_rdown,l) * de(ij+u_rdown) ) hv = Riv2(ij,vdown) * rhodz(ij,l) & + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) & + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv ENDDO ENDDO ENDDO !!! Compute potential vorticity contribution to du DO l=1,llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i du(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) * & ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+ & wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+ & wee(ij+u_right,3,1)*Fe(ij+u_left,l)+ & wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+ & wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+ & wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+ & wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+ & wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+ & wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+ & wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) ) du(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) * & ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+ & wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+ & wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+ & wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+ & wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+ & wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+ & wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+ & wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+ & wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+ & wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) ) du(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) * & ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+ & wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+ & wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+ & wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+ & wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+ & wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+ & wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+ & wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+ & wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+ & wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) ) ENDDO ENDDO ENDDO !!! Compute bernouilli term = Kinetic Energy + geopotential DO l=1,llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i berni(ij,l) = phi(ij,l) & + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 + & le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) ENDDO ENDDO ENDDO !!! second contribution to du DO l=1,llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * ( & 0.5*(theta(ij,l)+theta(ij+t_right,l)) & *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) & + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) ) du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * ( & 0.5*(theta(ij,l)+theta(ij+t_lup,l)) & *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l)) & + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) ) du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * ( & 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l)) & + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) ) ENDDO ENDDO ENDDO !!! second contribution to du DO l=1,llm !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i out_u(ij+u_right,l)= 1/de(ij+u_right) * ( & 0.5*(theta(ij,l)+theta(ij+t_right,l)) & *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) & + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) ) out_u(ij+u_lup,l)= 1/de(ij+u_lup) * ( & 0.5*(theta(ij,l)+theta(ij+t_lup,l)) & *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l)) & + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) ) out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * ( & 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l)) & + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) ) ENDDO ENDDO ENDDO !!! contribution due to vertical advection ! Contribution to dtheta DO l=1,llm-1 !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ww = 0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) ) dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - ww dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1) + ww ENDDO ENDDO ENDDO ! Contribution to du DO l=1,llm-1 !$OMP DO DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ww = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1)) uu = u(ij+u_right,l+1) - u(ij+u_right,l) du(ij+u_right, l ) = du(ij+u_right,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))) du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1))) ww = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1)) uu = u(ij+u_lup,l+1) - u(ij+u_lup,l) du(ij+u_lup, l ) = du(ij+u_lup,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))) du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1))) ww = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1)) uu = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) du(ij+u_ldown, l ) = du(ij+u_ldown,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))) du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1))) ENDDO ENDDO ENDDO !!$OMP BARRIER !!$OMP MASTER ! DEALLOCATE(theta) ! potential temperature ! DEALLOCATE(p) ! pression ! DEALLOCATE(pk) ! Exner function ! DEALLOCATE(pks) ! DEALLOCATE(alpha) ! DEALLOCATE(beta) ! DEALLOCATE(phi) ! geopotential ! DEALLOCATE(mass) ! mass ! DEALLOCATE(rhodz) ! mass density ! DEALLOCATE(Fe) ! mass flux ! DEALLOCATE(Ftheta) ! theta flux ! DEALLOCATE(convm) ! mass flux convergence ! DEALLOCATE(w) ! vertical velocity ! DEALLOCATE(qv) ! potential velocity ! DEALLOCATE(berni) ! bernouilli term !!$OMP END MASTER !!$OMP BARRIER END SUBROUTINE compute_caldyn END MODULE caldyn_gcm_mod