[906] | 1 | MODULE compute_rhodz_mod |
---|
| 2 | USE earth_const, ONLY : g |
---|
| 3 | USE disvert_mod, ONLY : ap, bp |
---|
| 4 | USE grid_param, ONLY : llm |
---|
| 5 | USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT |
---|
| 6 | IMPLICIT NONE |
---|
| 7 | PRIVATE |
---|
| 8 | SAVE |
---|
| 9 | |
---|
| 10 | #include "../unstructured/unstructured.h90" |
---|
| 11 | |
---|
| 12 | PUBLIC :: compute_rhodz_hex, compute_rhodz_unst |
---|
| 13 | |
---|
| 14 | CONTAINS |
---|
| 15 | |
---|
| 16 | #if BEGIN_DYSL |
---|
| 17 | |
---|
| 18 | KERNEL(compute_rhodz) |
---|
| 19 | IF(comp) THEN |
---|
| 20 | FORALL_CELLS() |
---|
| 21 | ON_PRIMAL |
---|
| 22 | m = ( AP(CELL)-AP(UP(CELL)) + (BP(CELL)-BP(UP(CELL)))*ps(HIDX(CELL)) )/g |
---|
| 23 | rhodz(CELL)=m |
---|
| 24 | END_BLOCK |
---|
| 25 | END_BLOCK |
---|
| 26 | ELSE |
---|
| 27 | err=0. |
---|
| 28 | FORALL_CELLS_EXT() |
---|
| 29 | ON_PRIMAL |
---|
| 30 | m = ( AP(CELL)-AP(UP(CELL)) + (BP(CELL)-BP(UP(CELL)))*ps(HIDX(CELL)) )/g |
---|
| 31 | err = MAX(err, ABS(m-rhodz(CELL))) |
---|
| 32 | END_BLOCK |
---|
| 33 | END_BLOCK |
---|
| 34 | IF(err>1e-10) THEN |
---|
| 35 | PRINT *, 'Discrepancy between ps and rhodz detected', err |
---|
| 36 | STOP |
---|
| 37 | END IF |
---|
| 38 | END IF |
---|
| 39 | END_BLOCK |
---|
| 40 | |
---|
| 41 | #endif END_DYSL |
---|
| 42 | |
---|
| 43 | SUBROUTINE compute_rhodz_unst(comp, ps, rhodz) |
---|
| 44 | USE data_unstructured_mod, ONLY : primal_num |
---|
| 45 | LOGICAL, INTENT(IN) :: comp |
---|
| 46 | FIELD_PS, INTENT(IN) :: ps |
---|
| 47 | FIELD_MASS, INTENT(INOUT) :: rhodz |
---|
| 48 | DECLARE_INDICES |
---|
| 49 | NUM :: m, err |
---|
| 50 | #define AP(l,ij) ap(l) |
---|
| 51 | #define BP(l,ij) bp(l) |
---|
| 52 | #include "../kernels_unst/compute_rhodz.k90" |
---|
| 53 | #undef AP |
---|
| 54 | #undef BP |
---|
| 55 | END SUBROUTINE compute_rhodz_unst |
---|
| 56 | |
---|
| 57 | SUBROUTINE compute_rhodz_hex(comp, ps, rhodz) |
---|
| 58 | USE icosa |
---|
| 59 | USE omp_para |
---|
| 60 | LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check |
---|
| 61 | REAL(rstd), INTENT(IN) :: ps(iim*jjm) |
---|
| 62 | REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) |
---|
| 63 | REAL(rstd) :: m, err |
---|
| 64 | INTEGER :: ij,l |
---|
| 65 | #define AP(ij,l) ap(l) |
---|
| 66 | #define BP(ij,l) bp(l) |
---|
| 67 | #include "../kernels_hex/compute_rhodz.k90" |
---|
| 68 | #undef AP |
---|
| 69 | #undef BP |
---|
| 70 | END SUBROUTINE compute_rhodz_hex |
---|
| 71 | |
---|
| 72 | SUBROUTINE compute_rhodz_handmade(comp, ps, rhodz) |
---|
| 73 | USE icosa |
---|
| 74 | USE omp_para |
---|
| 75 | LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check |
---|
| 76 | REAL(rstd), INTENT(IN) :: ps(iim*jjm) |
---|
| 77 | REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) |
---|
| 78 | REAL(rstd) :: m, err |
---|
| 79 | INTEGER :: l,i,j,ij,dd |
---|
| 80 | err=0. |
---|
| 81 | |
---|
| 82 | IF(comp) THEN |
---|
| 83 | dd=1 |
---|
| 84 | ELSE |
---|
| 85 | dd=0 |
---|
| 86 | END IF |
---|
| 87 | |
---|
| 88 | DO l = ll_begin, ll_end |
---|
| 89 | DO j=jj_begin-dd,jj_end+dd |
---|
| 90 | DO i=ii_begin-dd,ii_end+dd |
---|
| 91 | ij=(j-1)*iim+i |
---|
| 92 | m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g |
---|
| 93 | IF(comp) THEN |
---|
| 94 | rhodz(ij,l) = m |
---|
| 95 | ELSE |
---|
| 96 | err = MAX(err,abs(m-rhodz(ij,l))) |
---|
| 97 | END IF |
---|
| 98 | ENDDO |
---|
| 99 | ENDDO |
---|
| 100 | ENDDO |
---|
| 101 | |
---|
| 102 | IF(.NOT. comp) THEN |
---|
| 103 | IF(err>1e-10) THEN |
---|
| 104 | PRINT *, 'Discrepancy between ps and rhodz detected', err |
---|
| 105 | STOP |
---|
| 106 | ELSE |
---|
| 107 | ! PRINT *, 'No discrepancy between ps and rhodz detected' |
---|
| 108 | END IF |
---|
| 109 | END IF |
---|
| 110 | |
---|
| 111 | END SUBROUTINE compute_rhodz_handmade |
---|
| 112 | |
---|
| 113 | END MODULE compute_rhodz_mod |
---|