Changeset 878 for codes/icosagcm/devel/Python

Timestamp:

05/29/19 20:33:00 (5 years ago)

Author:

jisesh

Message:

devel: moved DYSL into compute_caldyn_solver.F90,compute_theta.F90 and compute_NH_geopot.F90

Location:

codes/icosagcm/devel/Python/src

Files:

• kernels_caldyn_NH.jin (modified) (2 diffs)
• kernels_caldyn_hevi.jin (modified) (1 diff)

codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin

@@ -10 +10 @@
 {%- endmacro %}
 
-KERNEL(compute_NH_geopot)
-  tau2_g=tau*tau/g
-  g2=g*g
-  gm2 = 1./g2
-  vreff = Treff*cpp/preff*kappa
-  gamma = 1./(1.-kappa)
-    
-  BARRIER
-
-  ! compute Phi_star
-  SEQUENCE_C1
-    BODY('1,llm+1')
-      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)
-    END_BLOCK
-  END_BLOCK
-
-  ! Newton-Raphson iteration : Phi_il contains current guess value
-  DO iter=1,2 ! 2 iterations should be enough
-    ! Compute pressure, A_ik
-    SELECT CASE(caldyn_thermo)
-    CASE(thermo_theta)
-      {% call() compute_p_and_Aik() %}
-        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij
-        p_ik(CELL) = preff*(X_ij**gamma)
-        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
-      {% endcall %}
-    CASE(thermo_entropy)
-      {% call() compute_p_and_Aik() %}
-        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp
-        p_ik(CELL) = preff*exp(X_ij*gamma)
-        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
-      {% endcall %}
-    CASE DEFAULT
-        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
-        STOP
-    END SELECT
-    
-    ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
-    ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
-
-    SEQUENCE_C1
-      ! Compute residual R_il and B_il
-      PROLOGUE(1)
-        ! bottom interface l=1
-        ml_g2 = gm2*m_il(CELL) 
-        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot
-        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
-                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )
-      END_BLOCK
-      BODY('2,llm')
-        ! inner interfaces
-        ml_g2 = gm2*m_il(CELL) 
-        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2
-        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
-                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))
-        ! consistency check : if Wil=0 and initial state is in hydrostatic balance
-        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2
-        ! and residual = tau^2*(ml+(1/g)dl_pi)=0
-      END_BLOCK
-      EPILOGUE('llm+1')
-        ! top interface l=llm+1
-        ml_g2 = gm2*m_il(CELL) 
-        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2
-        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
-                  + tau2_g*( ptop-p_ik(DOWN(CELL)) )
-      END_BLOCK
-      !
-      ! Forward sweep :
-      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 
-      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
-      PROLOGUE(1)
-        X_ij = 1./B_il(CELL)
-        C_ik(CELL) = -A_ik(CELL) * X_ij 
-        D_il(CELL) =  R_il(CELL) * X_ij
-      END_BLOCK
-      BODY('2,llm')
-        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
-        C_ik(CELL) = -A_ik(CELL) * X_ij 
-        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
-      END_BLOCK
-      EPILOGUE('llm+1')
-        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
-        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
-        ! Back substitution :
-        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
-        ! + Newton-Raphson update
-        ! top interface l=llm+1
-        x_il(CELL) = D_il(CELL)
-        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
-      END_BLOCK
-      BODY('llm,1,-1')
-        ! Back substitution at lower interfaces
-        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))
-        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
-      END_BLOCK
-    END_BLOCK
-
-    IF(debug_hevi_solver) THEN
-       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
-       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
-       DO l=1,llm+1
-          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
-       END DO
-       DO l=1,llm+1
-          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:))) 
-       END DO
-    END IF
-
-  END DO ! Newton-Raphson
-
-  BARRIER
-
-  debug_hevi_solver=.FALSE.
-END_BLOCK
-
 KERNEL(caldyn_mil)
   FORALL_CELLS_EXT('1', 'llm+1')
@@ -131 +16 @@
        CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
        CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
-    END_BLOCK
-  END_BLOCK
-END_BLOCK
-
-KERNEL(caldyn_solver)
-  ! 
-  ! Compute pressure (pres) and Exner function (pk)
-  ! kappa = R/Cp
-  ! 1-kappa = Cv/Cp
-  ! Cp/Cv = 1/(1-kappa)
-  gamma    = 1./(1.-kappa)
-  vreff    = Rd*Treff/preff ! reference specific volume
-  Cvd      = 1./(cpp-Rd)
-  Rd_preff = kappa*cpp/preff
-  FORALL_CELLS_EXT()
-    SELECT CASE(caldyn_thermo)
-    CASE(thermo_theta)
-      ON_PRIMAL
-        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
-        rho_ij = rho_ij*g*rhodz(CELL)
-        X_ij = Rd_preff*theta(CELL,1)*rho_ij
-        ! kappa.theta.rho = p/exner
-        ! => X = (p/p0)/(exner/Cp)
-        !      = (p/p0)^(1-kappa)
-        pres(CELL) = preff*(X_ij**gamma)  ! pressure
-        ! Compute Exner function (needed by compute_caldyn_fast)
-        ! other formulae possible if exponentiation is slow
-        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner
-      END_BLOCK
-    CASE(thermo_entropy)
-      ON_PRIMAL
-        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
-        rho_ij = rho_ij*g*rhodz(CELL)
-        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd )
-        pres(CELL) = rho_ij*Rd*T_ij
-        pk(CELL)   = T_ij
-      END_BLOCK
-    CASE DEFAULT
-      STOP
-    END SELECT
-  END_BLOCK
-
-  ! We need a barrier here because we compute pres above and do a vertical difference below
-  BARRIER
-
-  FORALL_CELLS_EXT('1', 'llm+1')
-    ON_PRIMAL
-      CST_IFTHEN(IS_BOTTOM_LEVEL)
-        ! Lower BC
-        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)
-      CST_ELSEIF(IS_TOP_INTERFACE)
-        ! Top BC
-        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)
-      CST_ELSE
-        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)
-      CST_ENDIF
-      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W
-      dPhi(CELL) = g*g*W(CELL)/m_il(CELL)
-    END_BLOCK
-  END_BLOCK
-
-  ! We need a barrier here because we update W above and do a vertical average below
-  BARRIER
-
-  FORALL_CELLS_EXT()
-    ON_PRIMAL
-      ! compute du = -0.5*g^2.grad(w^2)
-      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )
-    END_BLOCK
-  END_BLOCK
-  FORALL_CELLS_EXT()
-    ON_EDGES
-      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
     END_BLOCK
   END_BLOCK

codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin

@@ -96 +96 @@
 END_BLOCK
 
-KERNEL(theta)
-  IF(caldyn_eta==eta_mass) THEN ! Compute mass
-    ! FIXME : here mass_col is computed from rhodz
-    ! so that the DOFs are the same whatever caldyn_eta
-    ! in DYNAMICO mass_col is prognosed rather than rhodz
-    SEQUENCE_C1
-      PROLOGUE(0)
-        mass_col(HIDX(CELL))=0.
-      END_BLOCK
-      BODY('1,llm')
-          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
-      END_BLOCK
-    END_BLOCK
-    FORALL_CELLS_EXT()
-      ON_PRIMAL
-        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
-        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
-        !        rhodz(CELL) = m/g
-        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
-      END_BLOCK
-    END_BLOCK
-  END IF
-  DO iq=1,nqdyn
-    FORALL_CELLS_EXT()
-      ON_PRIMAL
-        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
-      END_BLOCK
-    END_BLOCK
-  END DO
-END_BLOCK