Changeset 937 for codes/icosagcm/devel/src
- Timestamp:
- 07/03/19 18:27:33 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/compute_geopot.F90
r920 r937 10 10 SAVE 11 11 12 PUBLIC :: compute_geopot_hex, compute_geopot_manual 12 #include "../unstructured/unstructured.h90" 13 14 PUBLIC :: compute_geopot_unst, compute_geopot_hex, compute_geopot_manual 13 15 14 16 CONTAINS 15 17 18 #ifdef BEGIN_DYSL 19 20 #define THECELL {{ thecell }} 21 22 {# ---------------- macro to generate code computing pressure top-down --------------- 23 formula = formula to compute 'gravitational' mass 24 = rhodz (dry) rhodz*theta (boussinesq) rhodz*(1+qv) (moist) #} 25 26 #define BALANCE(formula) {% call(thecell) balance() %} formula {% endcall %} 27 {% macro balance() %} {% set formula=caller %} 28 SEQUENCE_C1 29 PROLOGUE(llm) 30 pk(CELL) = ptop + .5*g*{{ formula('CELL') }} 31 END_BLOCK 32 BODY('llm-1,1,-1') 33 pk(CELL) = pk(UP(CELL)) + (.5*g)*({{ formula('CELL') }}+{{ formula('UP(CELL)') }}) 34 END_BLOCK 35 IF(caldyn_eta == eta_lag) THEN 36 EPILOGUE(1) 37 ps(HIDX(CELL)) = pk(CELL) + .5*g*{{ formula('CELL') }} 38 END_BLOCK 39 END IF 40 END_BLOCK 41 {%- endmacro %} 42 43 {# ------------ macro to generate code computing geopotential bottom-up -------------- 44 var = variable to be stored in pk(CELL) 45 caller() computes gv = g*v where v = specific volume 46 details depend on caldyn_thermo #} 47 48 #define GEOPOT(var) {% call geopot(var) %} 49 {% macro geopot(var) %} {% set formula=caller %} 50 SEQUENCE_C1 51 BODY('1,llm') 52 p_ik = pk(CELL) 53 {{ formula() }} 54 pk(CELL) = {{ var }} 55 geopot(UP(CELL)) = geopot(CELL) + gv*rhodz(CELL) 56 END_BLOCK 57 END_BLOCK 58 {%- endmacro %} 59 60 #define END_GEOPOT {% endcall %} 61 62 KERNEL(compute_geopot) 63 SELECT CASE(caldyn_thermo) 64 CASE(thermo_boussinesq) 65 ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure) 66 BALANCE( theta(THECELL,1)*rhodz(THECELL) ) 67 ! now pk contains the Lagrange multiplier (pressure) 68 ! specific volume 1 = dphi/g/rhodz 69 SEQUENCE_C1 70 BODY('1,llm') 71 geopot(UP(CELL)) = geopot(CELL) + g*rhodz(CELL) 72 END_BLOCK 73 END_BLOCK 74 CASE(thermo_theta) 75 BALANCE( rhodz(THECELL) ) 76 GEOPOT('exner_ik') 77 exner_ik = cpp * (p_ik/preff) ** kappa 78 gv = (g*kappa)*theta(CELL,1)*exner_ik/p_ik 79 END_GEOPOT 80 CASE(thermo_entropy) 81 BALANCE( rhodz(THECELL) ) 82 GEOPOT('temp_ik') 83 temp_ik = Treff*exp((theta(CELL,1) + Rd*log(p_ik/preff))/cpp) 84 gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p 85 END_GEOPOT 86 CASE(thermo_variable_Cp) 87 ! thermodynamics with variable Cp 88 ! Cp.dT = dh = Tds + vdp 89 ! pv = RT 90 ! => ds = (dh+v.dp)/T = Cp.dT/T - R dp/p 91 ! Cp(T) = Cp0 * (T/T0)^nu 92 ! => s(p,T) = Cp(T)/nu - R log(p/preff) 93 ! h = Cp(T).T/(nu+1) 94 BALANCE( rhodz(THECELL) ) 95 GEOPOT('temp_ik') 96 Cp_ik = nu*( theta(CELL,1) + Rd*log(p_ik/preff) ) 97 temp_ik = Treff* (Cp_ik/cpp)**(1./nu) 98 gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p 99 END_GEOPOT 100 CASE(thermo_moist) 101 BALANCE( rhodz(THECELL)*(1.+theta(THECELL,2)) ) 102 GEOPOT('temp_ik') 103 qv = theta(CELL,2) ! water vaper mixing ratio = mv/md 104 Rmix = Rd+qv*Rv 105 chi = ( theta(CELL,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) 106 temp_ik = Treff*exp(chi) 107 ! specific volume v = R*T/p 108 ! R = (Rd + qv.Rv)/(1+qv) 109 gv = g*Rmix*temp_ik/(p_ik*(1+qv)) 110 END_GEOPOT 111 END SELECT 112 END_BLOCK 113 114 #endif END_DYSL 115 16 116 !**************************** Geopotential ***************************** 17 117 118 SUBROUTINE compute_geopot_unst(rhodz,theta,ps,pk,geopot) 119 USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 120 USE data_unstructured_mod, ONLY : enter_trace, exit_trace, & 121 id_geopot 122 123 FIELD_MASS :: rhodz,pk ! IN, OUT 124 FIELD_THETA :: theta ! IN 125 FIELD_GEOPOT :: geopot ! IN(l=1)/OUT(l>1) 126 FIELD_PS :: ps ! OUT 127 DECLARE_INDICES 128 NUM :: chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix, Cp_ik 129 START_TRACE(id_geopot, 3,0,3) 130 #include "../kernels_unst/compute_geopot.k90" 131 STOP_TRACE 132 END SUBROUTINE compute_geopot_unst 133 18 134 SUBROUTINE compute_geopot_hex(rhodz,theta, ps,pk,geopot) 19 135 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
Note: See TracChangeset
for help on using the changeset viewer.