diff --git a/src/LIB/BITREP/modi_bitrep.f90 b/src/LIB/BITREP/modi_bitrep.f90 index de2471e0895e34e755c1c136e48691c06b81b6c2..783cca54db6e44737c2bad09bf7c415af25ed686 100644 --- a/src/LIB/BITREP/modi_bitrep.f90 +++ b/src/LIB/BITREP/modi_bitrep.f90 @@ -56,5 +56,15 @@ BR_POW = BR_EXP( PPOW * BR_LOG(PVAL) ) ! END FUNCTION ! +ELEMENTAL FUNCTION BR_P2(PVAL) +!$acc routine seq +! +REAL, INTENT(IN) :: PVAL +REAL :: BR_P2 +! +BR_P2 = PVAL * PVAL +!!$BR_P2 = PVAL ** 2 +! +END FUNCTION BR_P2 ! END MODULE MODI_BITREP diff --git a/src/MNH/advec_weno_k_2_aux.f90 b/src/MNH/advec_weno_k_2_aux.f90 index 3977e4e88073cb3aea1fffc251918bb46cf07103..ecd24a10254c5f91c8b6003ab1daaf08bb148ce3 100644 --- a/src/MNH/advec_weno_k_2_aux.f90 +++ b/src/MNH/advec_weno_k_2_aux.f90 @@ -298,6 +298,11 @@ USE MODD_LUNIT USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODI_GET_HALO +USE MODE_MPPDB +! +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -434,6 +439,7 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_UX::CYCL not yet tested' ZFNEG1(IE, :,:) = 0.5 * (3.0*PSRC(IE+1, :,:) - TEAST(:,:)) ZFNEG2(IW-1:IE,:,:) = 0.5 * (PSRC(IW-1:IE,:,:) + PSRC(IW:IE+1,:,:)) ZFNEG2(IE+1, :,:) = 0.5 * (PSRC(IE+1, :,:) + TEAST(:,:)) +#ifndef MNH_BITREP ! ! smoothness indicators for positive wind case ! @@ -444,7 +450,7 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_UX::CYCL not yet tested' ZBPOS2(IE+1, :,:) = (TEAST(:,:) - PSRC(IE+1, :,:))**2 ! ! smoothness indicators for negative wind case -! +! ZBNEG1(IW-1:IE-1,:,:) = (PSRC(IW:IE,:,:) - PSRC(IW+1:IE+1,:,:))**2 ZBNEG1(IE, :,:) = (PSRC(IE+1, :,:) - TEAST(:,:))**2 ZBNEG2(IW-1:IE,:,:) = (PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))**2 @@ -456,6 +462,30 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_UX::CYCL not yet tested' ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2 +#else +! +! smoothness indicators for positive wind case +! + ZBPOS1(IW:IE+1,:,:) = BR_P2(PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:)) + ZBPOS1(IW-1, :,:) = BR_P2(PSRC(IW-1, :,:) - TWEST(:,:)) +! + ZBPOS2(IW-1:IE,:,:) = BR_P2(PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:)) + ZBPOS2(IE+1, :,:) = BR_P2(TEAST(:,:) - PSRC(IE+1, :,:)) +! +! smoothness indicators for negative wind case +! + ZBNEG1(IW-1:IE-1,:,:) = BR_P2(PSRC(IW:IE,:,:) - PSRC(IW+1:IE+1,:,:)) + ZBNEG1(IE, :,:) = BR_P2(PSRC(IE+1, :,:) - TEAST(:,:)) + ZBNEG2(IW-1:IE,:,:) = BR_P2(PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:)) + ZBNEG2(IE+1, :,:) = BR_P2(PSRC(IE+1, :,:) - TEAST(:,:)) +! +! WENO weights +! + ZOMP1 = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1) + ZOMP2 = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2) + ZOMN1 = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1) + ZOMN2 = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2) +#endif ! ! WENO fluxes ! @@ -479,6 +509,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(IW-1,:,:) = 0.5 * (3.0*PSRC(IW-1,:,:) - TWEST(:,:)) ZFPOS2(IW-1,:,:) = 0.5 * (PSRC(IW-1, :,:) + PSRC(IW,:,:)) ZBPOS1(IW-1,:,:) = (PSRC(IW-1,:,:) - TWEST(:,:))**2 @@ -493,6 +524,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IW-1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW-1,:,:))**2 ZOMN1(IW-1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW-1,:,:))**2 ZOMN2(IW-1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW-1,:,:))**2 +#else + ZFPOS1(IW-1,:,:) = 0.5 * (3.0*PSRC(IW-1,:,:) - TWEST(:,:)) + ZFPOS2(IW-1,:,:) = 0.5 * (PSRC(IW-1, :,:) + PSRC(IW,:,:)) + ZBPOS1(IW-1,:,:) = BR_P2(PSRC(IW-1,:,:) - TWEST(:,:)) + ZBPOS2(IW-1,:,:) = BR_P2(PSRC(IW, :,:) - PSRC(IW-1,:,:)) +! + ZFNEG1(IW-1,:,:) = 0.5 * (3.0*PSRC(IW,:,:) - PSRC(IW+1,:,:)) + ZFNEG2(IW-1,:,:) = 0.5 * (PSRC(IW-1, :,:) + PSRC(IW, :,:)) + ZBNEG1(IW-1,:,:) = BR_P2(PSRC(IW, :,:) - PSRC(IW+1,:,:)) + ZBNEG2(IW-1,:,:) = BR_P2(PSRC(IW-1,:,:) - PSRC(IW, :,:)) +! + ZOMP1(IW-1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IW-1,:,:)) + ZOMP2(IW-1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IW-1,:,:)) + ZOMN1(IW-1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IW-1,:,:)) + ZOMN2(IW-1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IW-1,:,:)) +#endif ! PR(IW-1,:,:) = (ZOMN2(IW-1,:,:)/(ZOMN1(IW-1,:,:)+ZOMN2(IW-1,:,:)) * ZFNEG2(IW-1,:,:) + & (ZOMN1(IW-1,:,:)/(ZOMN1(IW-1,:,:)+ZOMN2(IW-1,:,:)) * ZFNEG1(IW-1,:,:))) * (0.5-SIGN(0.5,PRUCT(IW-1,:,:))) + & @@ -506,8 +553,9 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(IE,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:)) - ZFPOS2(IE,:,:) = 0.5 * (PSRC(IE, :,:) + PSRC(IE+1,:,:)) + ZFPOS2(IE,:,:) = 0.5 * (PSRC(IE, :,:) + PSRC(IE+1,:,:)) ZBPOS1(IE,:,:) = (PSRC(IE,:,:) - PSRC(IE-1,:,:))**2 ZBPOS2(IE,:,:) = (PSRC(IE+1,:,:) - PSRC(IE,:,:))**2 ! @@ -520,6 +568,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IE,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IE,:,:))**2 ZOMN1(IE,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IE,:,:))**2 ZOMN2(IE,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IE,:,:))**2 +#else + ZFPOS1(IE,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:)) + ZFPOS2(IE,:,:) = 0.5 * (PSRC(IE, :,:) + PSRC(IE+1,:,:)) + ZBPOS1(IE,:,:) = BR_P2(PSRC(IE,:,:) - PSRC(IE-1,:,:)) + ZBPOS2(IE,:,:) = BR_P2(PSRC(IE+1,:,:) - PSRC(IE,:,:)) +! + ZFNEG1(IE,:,:) = 0.5 * (3.0*PSRC(IE+1,:,:) - TEAST(:,:)) + ZFNEG2(IE,:,:) = 0.5 * (PSRC(IE,:,:) + PSRC(IE+1,:,:)) + ZBNEG1(IE,:,:) = BR_P2(PSRC(IE+1,:,:) - TEAST(:,:)) + ZBNEG2(IE,:,:) = BR_P2(PSRC(IE, :,:) - PSRC(IE+1,:,:)) +! + ZOMP1(IE,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IE,:,:)) + ZOMP2(IE,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IE,:,:)) + ZOMN1(IE,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IE,:,:)) + ZOMN2(IE,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IE,:,:)) +#endif ! PR(IE,:,:) = (ZOMN2(IE,:,:)/(ZOMN1(IE,:,:)+ZOMN2(IE,:,:)) * ZFNEG2(IE,:,:) + & (ZOMN1(IE,:,:)/(ZOMN1(IE,:,:)+ZOMN2(IE,:,:)) * ZFNEG1(IE,:,:))) * (0.5-SIGN(0.5,PRUCT(IE,:,:))) + & @@ -530,6 +594,7 @@ CASE ('OPEN','WALL','NEST') ! ! USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE ! +#ifndef MNH_BITREP ZFPOS1(IW:IE-1,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) ZFPOS2(IW:IE-1,:,:) = 0.5 * (PSRC(IW:IE-1, :,:) + PSRC(IW+1:IE, :,:)) ZBPOS1(IW:IE-1,:,:) = (PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))**2 @@ -544,6 +609,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IW:IE-1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW:IE-1,:,:))**2 ZOMN1(IW:IE-1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW:IE-1,:,:))**2 ZOMN2(IW:IE-1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW:IE-1,:,:))**2 +#else + ZFPOS1(IW:IE-1,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) + ZFPOS2(IW:IE-1,:,:) = 0.5 * (PSRC(IW:IE-1, :,:) + PSRC(IW+1:IE, :,:)) + ZBPOS1(IW:IE-1,:,:) = BR_P2(PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) + ZBPOS2(IW:IE-1,:,:) = BR_P2(PSRC(IW+1:IE,:,:) - PSRC(IW:IE-1, :,:)) +! + ZFNEG1(IW:IE-1,:,:) = 0.5 * (3.0*PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:)) + ZFNEG2(IW:IE-1,:,:) = 0.5 * (PSRC(IW:IE-1, :,:) + PSRC(IW+1:IE, :,:)) + ZBNEG1(IW:IE-1,:,:) = BR_P2(PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:)) + ZBNEG2(IW:IE-1,:,:) = BR_P2(PSRC(IW:IE-1,:,:) - PSRC(IW+1:IE,:,:)) +! + ZOMP1(IW:IE-1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IW:IE-1,:,:)) + ZOMP2(IW:IE-1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IW:IE-1,:,:)) + ZOMN1(IW:IE-1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IW:IE-1,:,:)) + ZOMN2(IW:IE-1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IW:IE-1,:,:)) +#endif ! PR(IW:IE-1,:,:) = (ZOMN2(IW:IE-1,:,:)/(ZOMN1(IW:IE-1,:,:)+ZOMN2(IW:IE-1,:,:)) * ZFNEG2(IW:IE-1,:,:) + & (ZOMN1(IW:IE-1,:,:)/(ZOMN1(IW:IE-1,:,:)+ZOMN2(IW:IE-1,:,:)) * ZFNEG1(IW:IE-1,:,:))) * (0.5-SIGN(0.5,PRUCT(IW:IE-1,:,:))) + & @@ -591,6 +672,9 @@ USE MODD_LUNIT USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -725,7 +809,8 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_MX::CYCL not yet tested' ! ZFNEG2(IW:IE+1,:,:) = 0.5 * (PSRC(IW:IE+1,:,:) + PSRC(IW-1:IE,:,:)) ZFNEG2(IW-1, :,:) = 0.5 * (PSRC(IW-1, :,:) + TWEST(:,:)) -! +#ifndef MNH_BITREP +! ! smoothness indicators for positive wind case ! ZBPOS1(IW+1:IE+1,:,:) = (PSRC(IW:IE,:,:) - PSRC(IW-1:IE-1,:,:))**2 @@ -736,7 +821,7 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_MX::CYCL not yet tested' ZBPOS2(IW-1, :,:) = (PSRC(IW-1, :,:) - TWEST(:,:))**2 ! ! smoothness indicators for negative wind case -! +! ZBNEG1(IW-1:IE,:,:) = (PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))**2 ZBNEG1(IE+1, :,:) = (PSRC(IE+1, :,:) - TEAST(:,:))**2 ! @@ -749,6 +834,32 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_MX::CYCL not yet tested' ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2 +#else +! +! smoothness indicators for positive wind case +! + ZBPOS1(IW+1:IE+1,:,:) = BR_P2(PSRC(IW:IE,:,:) - PSRC(IW-1:IE-1,:,:)) + ZBPOS1(IW, :,:) = BR_P2(PSRC(IW-1, :,:) - TWEST(:,:)) +!! ZBPOS1(IW-1, :,:) = - 999. +! + ZBPOS2(IW:IE+1,:,:) = BR_P2(PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:)) + ZBPOS2(IW-1, :,:) = BR_P2(PSRC(IW-1, :,:) - TWEST(:,:)) +! +! smoothness indicators for negative wind case +! + ZBNEG1(IW-1:IE,:,:) = BR_P2(PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:)) + ZBNEG1(IE+1, :,:) = BR_P2(PSRC(IE+1, :,:) - TEAST(:,:)) +! + ZBNEG2(IW:IE+1,:,:) = BR_P2(PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:)) + ZBNEG2(IW-1, :,:) = BR_P2(TWEST(:,:) - PSRC(IW-1,:,:)) +! +! WENO weights +! + ZOMP1 = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1) + ZOMP2 = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2) + ZOMN1 = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1) + ZOMN2 = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2) +#endif ! ! WENO fluxes ! @@ -772,6 +883,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(IW,:,:) = 0.5 * (3.0*PSRC(IW-1, :,:) - TWEST(:,:)) ZFPOS2(IW,:,:) = 0.5 * (PSRC(IW-1, :,:) + PSRC(IW, :,:)) ZBPOS1(IW,:,:) = (PSRC(IW-1,:,:) - TWEST(:,:))**2 @@ -786,6 +898,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IW,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW,:,:))**2 ZOMN1(IW,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW,:,:))**2 ZOMN2(IW,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW,:,:))**2 +#else + ZFPOS1(IW,:,:) = 0.5 * (3.0*PSRC(IW-1, :,:) - TWEST(:,:)) + ZFPOS2(IW,:,:) = 0.5 * (PSRC(IW-1, :,:) + PSRC(IW, :,:)) + ZBPOS1(IW,:,:) = BR_P2(PSRC(IW-1,:,:) - TWEST(:,:)) + ZBPOS2(IW,:,:) = BR_P2(PSRC(IW, :,:) - PSRC(IW-1,:,:)) +! + ZFNEG1(IW,:,:) = 0.5 * (3.0*PSRC(IW,:,:) - PSRC(IW+1,:,:)) + ZFNEG2(IW,:,:) = 0.5 * (PSRC(IW, :,:) + PSRC(IW-1,:,:)) + ZBNEG1(IW,:,:) = BR_P2(PSRC(IW,:,:) - PSRC(IW+1,:,:)) + ZBNEG2(IW,:,:) = BR_P2(PSRC(IW-1,:,:) - PSRC(IW,:,:)) +! + ZOMP1(IW,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IW,:,:)) + ZOMP2(IW,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IW,:,:)) + ZOMN1(IW,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IW,:,:)) + ZOMN2(IW,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IW,:,:)) +#endif ! PR(IW,:,:) = (ZOMP2(IW,:,:)/(ZOMP1(IW,:,:)+ZOMP2(IW,:,:)) * ZFPOS2(IW,:,:) + & (ZOMP1(IW,:,:)/(ZOMP1(IW,:,:)+ZOMP2(IW,:,:)) * ZFPOS1(IW,:,:))) * (0.5+SIGN(0.5,PRUCT(IW,:,:))) + & @@ -799,6 +927,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(IE+1,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:)) ZFPOS2(IE+1,:,:) = 0.5 * (PSRC(IE, :,:) + PSRC(IE+1,:,:)) ZBPOS1(IE+1,:,:) = (PSRC(IE,:,:) - PSRC(IE-1,:,:))**2 @@ -813,6 +942,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IE+1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IE+1,:,:))**2 ZOMN1(IE+1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IE+1,:,:))**2 ZOMN2(IE+1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IE+1,:,:))**2 +#else + ZFPOS1(IE+1,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:)) + ZFPOS2(IE+1,:,:) = 0.5 * (PSRC(IE, :,:) + PSRC(IE+1,:,:)) + ZBPOS1(IE+1,:,:) = BR_P2(PSRC(IE,:,:) - PSRC(IE-1,:,:)) + ZBPOS2(IE+1,:,:) = BR_P2(PSRC(IE+1,:,:) - PSRC(IE,:,:)) +! + ZFNEG1(IE+1,:,:) = 0.5 * (3.0*PSRC(IE+1,:,:) - TEAST(:,:)) + ZFNEG2(IE+1,:,:) = 0.5 * (PSRC(IE+1, :,:) + PSRC(IE,:,:)) + ZBNEG1(IE+1,:,:) = BR_P2(PSRC(IE+1,:,:) - TEAST(:,:)) + ZBNEG2(IE+1,:,:) = BR_P2(PSRC(IE, :,:) - PSRC(IE+1,:,:)) +! + ZOMP1(IE+1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IE+1,:,:)) + ZOMP2(IE+1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IE+1,:,:)) + ZOMN1(IE+1,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IE+1,:,:)) + ZOMN2(IE+1,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IE+1,:,:)) +#endif ! PR(IE+1,:,:) = (ZOMP2(IE+1,:,:)/(ZOMP1(IE+1,:,:)+ZOMP2(IE+1,:,:)) * ZFPOS2(IE+1,:,:) + & (ZOMP1(IE+1,:,:)/(ZOMP1(IE+1,:,:)+ZOMP2(IE+1,:,:)) * ZFPOS1(IE+1,:,:))) * (0.5+SIGN(0.5,PRUCT(IE+1,:,:))) + & @@ -823,6 +968,7 @@ CASE ('OPEN','WALL','NEST') ! ! USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE ! +#ifndef MNH_BITREP ZFPOS1(IW+1:IE,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) ZFPOS2(IW+1:IE,:,:) = 0.5 * (PSRC(IW:IE-1, :,:) + PSRC(IW+1:IE, :,:)) ZBPOS1(IW+1:IE,:,:) = (PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))**2 @@ -837,6 +983,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(IW+1:IE,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW+1:IE,:,:))**2 ZOMN1(IW+1:IE,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW+1:IE,:,:))**2 ZOMN2(IW+1:IE,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW+1:IE,:,:))**2 +#else + ZFPOS1(IW+1:IE,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) + ZFPOS2(IW+1:IE,:,:) = 0.5 * (PSRC(IW:IE-1, :,:) + PSRC(IW+1:IE, :,:)) + ZBPOS1(IW+1:IE,:,:) = BR_P2(PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:)) + ZBPOS2(IW+1:IE,:,:) = BR_P2(PSRC(IW+1:IE,:,:) - PSRC(IW:IE-1,:,:)) +! + ZFNEG1(IW+1:IE,:,:) = 0.5 * (3.0*PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:)) + ZFNEG2(IW+1:IE,:,:) = 0.5 * (PSRC(IW+1:IE, :,:) + PSRC(IW:IE-1, :,:)) + ZBNEG1(IW+1:IE,:,:) = BR_P2(PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:)) + ZBNEG2(IW+1:IE,:,:) = BR_P2(PSRC(IW:IE-1,:,:) - PSRC(IW+1:IE,:,:)) +! + ZOMP1(IW+1:IE,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(IW+1:IE,:,:)) + ZOMP2(IW+1:IE,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(IW+1:IE,:,:)) + ZOMN1(IW+1:IE,:,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(IW+1:IE,:,:)) + ZOMN2(IW+1:IE,:,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(IW+1:IE,:,:)) +#endif ! PR(IW+1:IE,:,:) = (ZOMP2(IW+1:IE,:,:)/(ZOMP1(IW+1:IE,:,:)+ZOMP2(IW+1:IE,:,:)) * ZFPOS2(IW+1:IE,:,:) + & (ZOMP1(IW+1:IE,:,:)/(ZOMP1(IW+1:IE,:,:)+ZOMP2(IW+1:IE,:,:)) * ZFPOS1(IW+1:IE,:,:))) * (0.5+SIGN(0.5,PRUCT(IW+1:IE,:,:))) + & @@ -884,6 +1046,9 @@ USE MODD_LUNIT USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -1017,12 +1182,13 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_MY::CYCL not yet tested' ! ZFNEG2(:,IS:IN+1,:) = 0.5 * (PSRC(:,IS:IN+1,:) + PSRC(:,IS-1:IN,:)) ZFNEG2(:,IS-1, :) = 0.5 * (PSRC(:,IS-1, :) + TSOUTH(:,:)) +#ifndef MNH_BITREP ! ! smoothness indicators for positive wind case ! ZBPOS1(:,IS+1:IN+1,:) = (PSRC(:,IS:IN,:) - PSRC(:,IS-1:IN-1,:))**2 ZBPOS1(:,IS, :) = (PSRC(:,IS-1, :) - TSOUTH(:,:))**2 -!! ZBPOS1(:,IS-1, :) = - 999. +!! ZBPOS1(:,IS-1, :) = - 999. ! ZBPOS2(:,IS:IN+1,:) = (PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:))**2 ZBPOS2(:,IS-1, :) = (PSRC(:,IS-1, :) - TSOUTH(:,:))**2 @@ -1041,6 +1207,32 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_MY::CYCL not yet tested' ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2 +#else +! +! smoothness indicators for positive wind case +! + ZBPOS1(:,IS+1:IN+1,:) = BR_P2(PSRC(:,IS:IN,:) - PSRC(:,IS-1:IN-1,:)) + ZBPOS1(:,IS, :) = BR_P2(PSRC(:,IS-1, :) - TSOUTH(:,:)) +!! ZBPOS1(:,IS-1, :) = - 999. +! + ZBPOS2(:,IS:IN+1,:) = BR_P2(PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:)) + ZBPOS2(:,IS-1, :) = BR_P2(PSRC(:,IS-1, :) - TSOUTH(:,:)) +! +! smoothness indicators for negative wind case +! + ZBNEG1(:,IS-1:IN,:) = BR_P2(PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:)) + ZBNEG1(:,IN+1, :) = BR_P2(PSRC(:,IN+1, :) - TNORTH(:,:)) +! + ZBNEG2(:,IS:IN+1,:) = BR_P2(PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:)) + ZBNEG2(:,IS-1, :) = BR_P2(TSOUTH(:,:) - PSRC(:,IS-1,:)) +! +! WENO weights +! + ZOMP1 = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1) + ZOMP2 = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2) + ZOMN1 = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1) + ZOMN2 = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2) +#endif ! ! WENO fluxes ! @@ -1064,6 +1256,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(:,IS,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TSOUTH(:,:)) ZFPOS2(:,IS,:) = 0.5 * (PSRC(:,IS-1,:) + PSRC(:,IS,:)) ZBPOS1(:,IS,:) = (PSRC(:,IS-1,:) - TSOUTH(:,:))**2 @@ -1078,6 +1271,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IS,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS,:))**2 ZOMN1(:,IS,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS,:))**2 ZOMN2(:,IS,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS,:))**2 +#else + ZFPOS1(:,IS,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TSOUTH(:,:)) + ZFPOS2(:,IS,:) = 0.5 * (PSRC(:,IS-1,:) + PSRC(:,IS,:)) + ZBPOS1(:,IS,:) = BR_P2(PSRC(:,IS-1,:) - TSOUTH(:,:)) + ZBPOS2(:,IS,:) = BR_P2(PSRC(:,IS, :) - PSRC(:,IS-1,:)) +! + ZFNEG1(:,IS,:) = 0.5 * (3.0*PSRC(:,IS,:) - PSRC(:,IS+1,:)) + ZFNEG2(:,IS,:) = 0.5 * (PSRC(:,IS, :) + PSRC(:,IS-1,:)) + ZBNEG1(:,IS,:) = BR_P2(PSRC(:,IS, :) - PSRC(:,IS+1,:)) + ZBNEG2(:,IS,:) = BR_P2(PSRC(:,IS-1,:) - PSRC(:,IS, :)) +! + ZOMP1(:,IS,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IS,:)) + ZOMP2(:,IS,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IS,:)) + ZOMN1(:,IS,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IS,:)) + ZOMN2(:,IS,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IS,:)) +#endif ! PR(:,IS,:) = (ZOMP2(:,IS,:)/(ZOMP1(:,IS,:)+ZOMP2(:,IS,:)) * ZFPOS2(:,IS,:) + & (ZOMP1(:,IS,:)/(ZOMP1(:,IS,:)+ZOMP2(:,IS,:)) * ZFPOS1(:,IS,:))) * (0.5+SIGN(0.5,PRVCT(:,IS,:))) + & @@ -1091,6 +1300,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(:,IN+1,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:)) ZFPOS2(:,IN+1,:) = 0.5 * (PSRC(:,IN, :) + PSRC(:,IN+1,:)) ZBPOS1(:,IN+1,:) = (PSRC(:,IN,:) - PSRC(:,IN-1,:))**2 @@ -1105,6 +1315,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IN+1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IN+1,:))**2 ZOMN1(:,IN+1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IN+1,:))**2 ZOMN2(:,IN+1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IN+1,:))**2 +#else + ZFPOS1(:,IN+1,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:)) + ZFPOS2(:,IN+1,:) = 0.5 * (PSRC(:,IN, :) + PSRC(:,IN+1,:)) + ZBPOS1(:,IN+1,:) = BR_P2(PSRC(:,IN,:) - PSRC(:,IN-1,:)) + ZBPOS2(:,IN+1,:) = BR_P2(PSRC(:,IN+1,:) - PSRC(:,IN,:)) +! + ZFNEG1(:,IN+1,:) = 0.5 * (3.0*PSRC(:,IN+1,:) - TNORTH(:,:)) + ZFNEG2(:,IN+1,:) = 0.5 * (PSRC(:,IN+1, :) + PSRC(:,IN,:)) + ZBNEG1(:,IN+1,:) = BR_P2(PSRC(:,IN+1,:) - TNORTH(:,:)) + ZBNEG2(:,IN+1,:) = BR_P2(PSRC(:,IN, :) - PSRC(:,IN+1,:)) +! + ZOMP1(:,IN+1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IN+1,:)) + ZOMP2(:,IN+1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IN+1,:)) + ZOMN1(:,IN+1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IN+1,:)) + ZOMN2(:,IN+1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IN+1,:)) +#endif ! PR(:,IN+1,:) = (ZOMP2(:,IN+1,:)/(ZOMP1(:,IN+1,:)+ZOMP2(:,IN+1,:)) * ZFPOS2(:,IN+1,:) + & (ZOMP1(:,IN+1,:)/(ZOMP1(:,IN+1,:)+ZOMP2(:,IN+1,:)) * ZFPOS1(:,IN+1,:))) * (0.5+SIGN(0.5,PRVCT(:,IN+1,:))) + & @@ -1115,6 +1341,7 @@ CASE ('OPEN','WALL','NEST') ! ! USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE ! +#ifndef MNH_BITREP ZFPOS1(:,IS+1:IN,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) ZFPOS2(:,IS+1:IN,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) ZBPOS1(:,IS+1:IN,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))**2 @@ -1129,6 +1356,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IS+1:IN,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS+1:IN,:))**2 ZOMN1(:,IS+1:IN,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS+1:IN,:))**2 ZOMN2(:,IS+1:IN,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS+1:IN,:))**2 +#else + ZFPOS1(:,IS+1:IN,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) + ZFPOS2(:,IS+1:IN,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) + ZBPOS1(:,IS+1:IN,:) = BR_P2(PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) + ZBPOS2(:,IS+1:IN,:) = BR_P2(PSRC(:,IS+1:IN,:) - PSRC(:,IS:IN-1, :)) +! + ZFNEG1(:,IS+1:IN,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) + ZFNEG2(:,IS+1:IN,:) = 0.5 * (PSRC(:,IS+1:IN, :) + PSRC(:,IS:IN-1, :)) + ZBNEG1(:,IS+1:IN,:) = BR_P2(PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) + ZBNEG2(:,IS+1:IN,:) = BR_P2(PSRC(:,IS:IN-1,:) - PSRC(:,IS+1:IN,:)) +! + ZOMP1(:,IS+1:IN,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IS+1:IN,:)) + ZOMP2(:,IS+1:IN,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IS+1:IN,:)) + ZOMN1(:,IS+1:IN,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IS+1:IN,:)) + ZOMN2(:,IS+1:IN,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IS+1:IN,:)) +#endif ! PR(:,IS+1:IN,:) = (ZOMP2(:,IS+1:IN,:)/(ZOMP1(:,IS+1:IN,:)+ZOMP2(:,IS+1:IN,:)) * ZFPOS2(:,IS+1:IN,:) + & (ZOMP1(:,IS+1:IN,:)/(ZOMP1(:,IS+1:IN,:)+ZOMP2(:,IS+1:IN,:)) * ZFPOS1(:,IS+1:IN,:))) * (0.5+SIGN(0.5,PRVCT(:,IS+1:IN,:))) + & @@ -1174,6 +1417,9 @@ USE MODD_LUNIT USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -1306,6 +1552,7 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_VY::CYCL not yet tested' ! ZFNEG2(:,IS-1:IN,:) = 0.5 * (PSRC(:,IS-1:IN,:) + PSRC(:,IS:IN+1,:)) ZFNEG2(:,IN+1, :) = 0.5 * (PSRC(:,IN+1, :) + TNORTH(:,:)) +#ifndef MNH_BITREP ! ! smoothness indicators for positive wind case ! @@ -1329,6 +1576,31 @@ PRINT *,'OPENACC: advec_weno_k_2_aux::ADVEC_WENO_K_2_VY::CYCL not yet tested' ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2 +#else +! +! smoothness indicators for positive wind case +! + ZBPOS1(:,IS:IN+1,:) = BR_P2(PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:)) + ZBPOS1(:,IS-1, :) = BR_P2(PSRC(:,IS-1, :) - TSOUTH(:,:)) +! + ZBPOS2(:,IS-1:IN,:) = BR_P2(PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:)) + ZBPOS2(:,IN+1, :) = BR_P2(TNORTH(:,:) - PSRC(:,IN+1, :)) +! +! smoothness indicators for negative wind case +! + ZBNEG1(:,IS-1:IN-1,:) = BR_P2(PSRC(:,IS:IN,:) - PSRC(:,IS+1:IN+1,:)) + ZBNEG1(:,IN, :) = BR_P2(PSRC(:,IN+1, :) - TNORTH(:,:)) +! + ZBNEG2(:,IS-1:IN,:) = BR_P2(PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:)) + ZBNEG2(:,IN+1, :) = BR_P2(PSRC(:,IN+1, :) - TNORTH(:,:)) +! +! WENO weights +! + ZOMP1 = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1) + ZOMP2 = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2) + ZOMN1 = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1) + ZOMN2 = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2) +#endif ! PR = (ZOMP2/(ZOMP1+ZOMP2) * ZFPOS2 + & (ZOMP1/(ZOMP1+ZOMP2) * ZFPOS1)) * (0.5+SIGN(0.5,PRVCT)) + & @@ -1350,6 +1622,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(:,IS-1,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TSOUTH(:,:)) ZFPOS2(:,IS-1,:) = 0.5 * (PSRC(:,IS-1, :) + PSRC(:,IS,:)) ZBPOS1(:,IS-1,:) = (PSRC(:,IS-1,:) - TSOUTH(:,:))**2 @@ -1364,6 +1637,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IS-1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS-1,:))**2 ZOMN1(:,IS-1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS-1,:))**2 ZOMN2(:,IS-1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS-1,:))**2 +#else + ZFPOS1(:,IS-1,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TSOUTH(:,:)) + ZFPOS2(:,IS-1,:) = 0.5 * (PSRC(:,IS-1, :) + PSRC(:,IS,:)) + ZBPOS1(:,IS-1,:) = BR_P2(PSRC(:,IS-1,:) - TSOUTH(:,:)) + ZBPOS2(:,IS-1,:) = BR_P2(PSRC(:,IS, :) - PSRC(:,IS-1,:)) +! + ZFNEG1(:,IS-1,:) = 0.5 * (3.0*PSRC(:,IS,:) - PSRC(:,IS+1,:)) + ZFNEG2(:,IS-1,:) = 0.5 * (PSRC(:,IS-1, :) + PSRC(:,IS,:)) + ZBNEG1(:,IS-1,:) = BR_P2(PSRC(:,IS,:) - PSRC(:,IS+1,:)) + ZBNEG2(:,IS-1,:) = BR_P2(PSRC(:,IS-1,:) - PSRC(:,IS,:)) +! + ZOMP1(:,IS-1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IS-1,:)) + ZOMP2(:,IS-1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IS-1,:)) + ZOMN1(:,IS-1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IS-1,:)) + ZOMN2(:,IS-1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IS-1,:)) +#endif ! PR(:,IS-1,:) = (ZOMP2(:,IS-1,:)/(ZOMP1(:,IS-1,:)+ZOMP2(:,IS-1,:)) * ZFPOS2(:,IS-1,:) + & (ZOMP1(:,IS-1,:)/(ZOMP1(:,IS-1,:)+ZOMP2(:,IS-1,:)) * ZFPOS1(:,IS-1,:))) * (0.5+SIGN(0.5,PRVCT(:,IS-1,:))) + & @@ -1377,6 +1666,7 @@ CASE ('OPEN','WALL','NEST') ! !!$ ELSEIF (NHALO == 1) THEN ELSE +#ifndef MNH_BITREP ZFPOS1(:,IN,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:)) ZFPOS2(:,IN,:) = 0.5 * (PSRC(:,IN, :) + PSRC(:,IN+1,:)) ZBPOS1(:,IN,:) = (PSRC(:,IN, :) - PSRC(:,IN-1,:))**2 @@ -1391,6 +1681,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IN,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IN,:))**2 ZOMN1(:,IN,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IN,:))**2 ZOMN2(:,IN,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IN,:))**2 +#else + ZFPOS1(:,IN,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:)) + ZFPOS2(:,IN,:) = 0.5 * (PSRC(:,IN, :) + PSRC(:,IN+1,:)) + ZBPOS1(:,IN,:) = BR_P2(PSRC(:,IN, :) - PSRC(:,IN-1,:)) + ZBPOS2(:,IN,:) = BR_P2(PSRC(:,IN+1,:) - PSRC(:,IN, :)) +! + ZFNEG1(:,IN,:) = 0.5 * (3.0*PSRC(:,IN+1,:) - TNORTH(:,:)) + ZFNEG2(:,IN,:) = 0.5 * (PSRC(:,IN, :) + PSRC(:,IN+1,:)) + ZBNEG1(:,IN,:) = BR_P2(PSRC(:,IN+1,:) - TNORTH(:,:)) + ZBNEG2(:,IN,:) = BR_P2(PSRC(:,IN, :) - PSRC(:,IN+1,:)) +! + ZOMP1(:,IN,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IN,:)) + ZOMP2(:,IN,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IN,:)) + ZOMN1(:,IN,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IN,:)) + ZOMN2(:,IN,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IN,:)) +#endif ! PR(:,IN,:) = (ZOMP2(:,IN,:)/(ZOMP1(:,IN,:)+ZOMP2(:,IN,:)) * ZFPOS2(:,IN,:) + & (ZOMP1(:,IN,:)/(ZOMP1(:,IN,:)+ZOMP2(:,IN,:)) * ZFPOS1(:,IN,:))) * (0.5+SIGN(0.5,PRVCT(:,IN,:))) + & @@ -1401,12 +1707,13 @@ CASE ('OPEN','WALL','NEST') ! ! USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE ! +#ifndef MNH_BITREP ZFPOS1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) ZFPOS2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) ZBPOS1(:,IS:IN-1,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))**2 ZBPOS2(:,IS:IN-1,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS:IN-1, :))**2 -! - ZFNEG1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) +! + ZFNEG1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) ZFNEG2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) ZBNEG1(:,IS:IN-1,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:))**2 ZBNEG2(:,IS:IN-1,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS+1:IN, :))**2 @@ -1415,6 +1722,22 @@ CASE ('OPEN','WALL','NEST') ZOMP2(:,IS:IN-1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS:IN-1,:))**2 ZOMN1(:,IS:IN-1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS:IN-1,:))**2 ZOMN2(:,IS:IN-1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS:IN-1,:))**2 +#else + ZFPOS1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) + ZFPOS2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) + ZBPOS1(:,IS:IN-1,:) = BR_P2(PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:)) + ZBPOS2(:,IS:IN-1,:) = BR_P2(PSRC(:,IS+1:IN,:) - PSRC(:,IS:IN-1, :)) +! + ZFNEG1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) + ZFNEG2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1, :) + PSRC(:,IS+1:IN, :)) + ZBNEG1(:,IS:IN-1,:) = BR_P2(PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:)) + ZBNEG2(:,IS:IN-1,:) = BR_P2(PSRC(:,IS:IN-1,:) - PSRC(:,IS+1:IN, :)) +! + ZOMP1(:,IS:IN-1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,IS:IN-1,:)) + ZOMP2(:,IS:IN-1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,IS:IN-1,:)) + ZOMN1(:,IS:IN-1,:) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,IS:IN-1,:)) + ZOMN2(:,IS:IN-1,:) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,IS:IN-1,:)) +#endif ! PR(:,IS:IN-1,:) = (ZOMP2(:,IS:IN-1,:)/(ZOMP1(:,IS:IN-1,:)+ZOMP2(:,IS:IN-1,:)) * ZFPOS2(:,IS:IN-1,:) + & (ZOMP1(:,IS:IN-1,:)/(ZOMP1(:,IS:IN-1,:)+ZOMP2(:,IS:IN-1,:)) * ZFPOS1(:,IS:IN-1,:))) * (0.5+SIGN(0.5,PRVCT(:,IS:IN-1,:))) + & @@ -1460,6 +1783,9 @@ USE MODE_ll USE MODD_CONF USE MODD_PARAMETERS,ONLY: JPVEXT USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -1518,10 +1844,11 @@ REAL, PARAMETER :: ZEPS = 1.0E-15 !* 0.3. COMPUTES THE DOMAIN DIMENSIONS ! ------------------------------ ! -!$acc kernels IB = 1 + JPVEXT IT = SIZE(PSRC,3) - JPVEXT ! +!$acc kernels +! PR(:,:,:) = 0.0 ! ZFPOS1 = 0.0 @@ -1550,6 +1877,7 @@ ZFPOS2(:,:,IB:IT-1) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT)) ! ZFNEG1(:,:,IB-1:IT-1) = 0.5 * (3.0*PSRC(:,:,IB:IT) - PSRC(:,:,IB+1:IT+1)) ZFNEG2(:,:,IB-1:IT) = 0.5 * (PSRC(:,:,IB-1:IT) + PSRC(:,:,IB:IT+1)) +#ifndef MNH_BITREP ! ! smoothness indicators for positive wind case ! @@ -1567,6 +1895,25 @@ ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2 ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2 +#else +! +! smoothness indicators for positive wind case +! +ZBPOS1(:,:,IB:IT-1) = BR_P2(PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2)) +ZBPOS2(:,:,IB:IT-1) = BR_P2(PSRC(:,:,IB+1:IT) - PSRC(:,:,IB:IT-1)) +! +! smoothness indicators for negative wind case +! +ZBNEG1(:,:,IB-1:IT-1) = BR_P2(PSRC(:,:,IB:IT) - PSRC(:,:,IB+1:IT+1)) +ZBNEG2(:,:,IB-1:IT) = BR_P2(PSRC(:,:,IB-1:IT) - PSRC(:,:,IB:IT+1)) +! +! WENO weights +! +ZOMP1 = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1) +ZOMP2 = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2) +ZOMN1 = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1) +ZOMN2 = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2) +#endif ! ! WENO fluxes ! @@ -1627,6 +1974,9 @@ USE MODE_ll USE MODD_CONF USE MODD_PARAMETERS,ONLY: JPVEXT USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -1687,10 +2037,11 @@ REAL, PARAMETER :: ZEPS = 1.0E-15 !* 0.3. COMPUTES THE DOMAIN DIMENSIONS ! ------------------------------ ! -!$acc kernels IB = 1 + JPVEXT IT = SIZE(PSRC,3) - JPVEXT ! +!$acc kernels +! PR(:,:,:) = 0.0 ! ZFPOS1 = 0.0 @@ -1717,6 +2068,7 @@ ZFPOS2(:,:,IB+1:IT) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT)) ! ZFNEG1(:,:,IB+1:IT) = 0.5 * (3.0*PSRC(:,:,IB+1:IT) - PSRC(:,:,IB+2:IT+1)) ZFNEG2(:,:,IB+1:IT) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT)) +#ifndef MNH_BITREP ! ! smoothness indicators for positive wind case ! @@ -1734,6 +2086,25 @@ ZOMP1(:,:,IB+1:IT) = ZGAMMA1 / (ZEPS + ZBPOS1(:,:,IB+1:IT))**2 ZOMP2(:,:,IB+1:IT) = ZGAMMA2 / (ZEPS + ZBPOS2(:,:,IB+1:IT))**2 ZOMN1(:,:,IB+1:IT) = ZGAMMA1 / (ZEPS + ZBNEG1(:,:,IB+1:IT))**2 ZOMN2(:,:,IB+1:IT) = ZGAMMA2 / (ZEPS + ZBNEG2(:,:,IB+1:IT))**2 +#else +! +! smoothness indicators for positive wind case +! +ZBPOS1(:,:,IB+1:IT) = BR_P2(PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2)) +ZBPOS2(:,:,IB+1:IT) = BR_P2(PSRC(:,:,IB+1:IT) - PSRC(:,:,IB:IT-1)) +! +! smoothness indicators for negative wind case +! +ZBNEG1(:,:,IB+1:IT) = BR_P2(PSRC(:,:,IB+1:IT) - PSRC(:,:,IB+2:IT+1)) +ZBNEG2(:,:,IB+1:IT) = BR_P2(PSRC(:,:,IB:IT-1) - PSRC(:,:,IB+1:IT)) +! +! WENO weights +! +ZOMP1(:,:,IB+1:IT) = ZGAMMA1 / BR_P2(ZEPS + ZBPOS1(:,:,IB+1:IT)) +ZOMP2(:,:,IB+1:IT) = ZGAMMA2 / BR_P2(ZEPS + ZBPOS2(:,:,IB+1:IT)) +ZOMN1(:,:,IB+1:IT) = ZGAMMA1 / BR_P2(ZEPS + ZBNEG1(:,:,IB+1:IT)) +ZOMN2(:,:,IB+1:IT) = ZGAMMA2 / BR_P2(ZEPS + ZBNEG2(:,:,IB+1:IT)) +#endif ! PR(:,:,IB+1:IT) = (ZOMP2(:,:,IB+1:IT)/(ZOMP1(:,:,IB+1:IT)+ZOMP2(:,:,IB+1:IT))* & ZFPOS2(:,:,IB+1:IT) + & diff --git a/src/MNH/advection_metsv.f90 b/src/MNH/advection_metsv.f90 index 5d409068496024e40f3fbf944d0801f470b6d9a2..1cca6070a69f67a277b33caa77d9cffd820a6993 100644 --- a/src/MNH/advection_metsv.f90 +++ b/src/MNH/advection_metsv.f90 @@ -178,6 +178,9 @@ USE MODE_FMWRIT USE MODE_DEVICE USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D #endif +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! !------------------------------------------------------------------------------- ! @@ -391,15 +394,27 @@ IF (.NOT. L1D) THEN ZCFLU(IIB:IIE,IJB:IJE,:) = ABS(ZRUCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP) ZCFLV(IIB:IIE,IJB:IJE,:) = ABS(ZRVCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP) ZCFLW(IIB:IIE,IJB:IJE,:) = ABS(ZRWCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP) +#ifndef MNH_BITREP IF (.NOT. L2D) THEN ZCFL(:,:,:) = SQRT(ZCFLU(:,:,:)**2+ZCFLV(:,:,:)**2+ZCFLW(:,:,:)**2) ELSE ZCFL = SQRT(ZCFLU(:,:,:)**2+ZCFLW(:,:,:)**2) END IF +#else + IF (.NOT. L2D) THEN + ZCFL(:,:,:) = SQRT(BR_P2(ZCFLU(:,:,:))+BR_P2(ZCFLV(:,:,:))+BR_P2(ZCFLW(:,:,:))) + ELSE + ZCFL = SQRT(BR_P2(ZCFLU(:,:,:))+BR_P2(ZCFLW(:,:,:))) + END IF +#endif ELSE ZCFLU = 0.0 ; ZCFLV = 0.0 ; ZCFLW = 0.0 ZCFLW(IIB:IIE,IJB:IJE,:) = ABS(ZRWCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP) +#ifndef MNH_BITREP ZCFL = SQRT(ZCFLW**2) +#else + ZCFL = SQRT(BR_P2(ZCFLW)) +#endif END IF !$acc end kernels ! diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90 index c710f8c34a89709b9a391c78bd7a0b879d1865ba..a36c445c6e07634880aa86a8187023d12eb07936 100644 --- a/src/MNH/bl89.f90 +++ b/src/MNH/bl89.f90 @@ -78,6 +78,9 @@ USE MODD_CONF, ONLY: CPROGRAM USE MODD_CST USE MODD_CTURB USE MODD_PARAMETERS +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! ! IMPLICIT NONE @@ -240,7 +243,11 @@ DO JK=IKTB,IKTE ZLWORK2= ( + ZG_O_THVREF(J1D,JK) * & ( ZVPT(J1D,JKK) - ZVPT(J1D,JK) ) & + SQRT (ABS( & +#ifndef MNH_BITREP ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - ZVPT(J1D,JK)) )**2 & +#else + BR_P2( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - ZVPT(J1D,JK)) ) & +#endif + 2. * ZINTE(J1D) * ZG_O_THVREF(J1D,JK) & * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK) ))) / & ( ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK)) @@ -278,7 +285,11 @@ DO JK=IKTB,IKTE ZLWORK2= ( - ZG_O_THVREF(J1D,JK) * & ( ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK) ) & + SQRT (ABS( & +#ifndef MNH_BITREP ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK)) )**2 & +#else + BR_P2( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK)) ) & +#endif + 2. * ZINTE(J1D) * ZG_O_THVREF(J1D,JK) & * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK) )) ) / & ( ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK) ) diff --git a/src/MNH/ground_paramn.f90 b/src/MNH/ground_paramn.f90 index 7c9ad622c52c5fa7779f3f354ea4c66f2482076e..c4b73c2c52d913d79ba5992da8ffb7a24601c667 100644 --- a/src/MNH/ground_paramn.f90 +++ b/src/MNH/ground_paramn.f90 @@ -164,6 +164,9 @@ USE MODD_TIME USE MODI_TEMPORAL_DIST ! USE MODD_PARAM_LIMA, ONLY : MSEDC=>LSEDC +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -397,7 +400,11 @@ CALL ROTATE_WIND(XUT,XVT,XWT, & ! 1.4 zonal and meridian components of the wind parallel to the slope ! --------------------------------------------------------------- ! +#ifndef MNH_BITREP ZWIND(:,:) = SQRT( ZUA**2 + ZVA**2 ) +#else +ZWIND(:,:) = SQRT( BR_P2(ZUA) + BR_P2(ZVA) ) +#endif ! ZU(:,:) = ZWIND(:,:) * SIN(ZDIR) ZV(:,:) = ZWIND(:,:) * COS(ZDIR) @@ -574,8 +581,13 @@ PSFU(:,:) = 0. PSFV(:,:) = 0. ! WHERE (ZSFU(:,:)/=XUNDEF .AND. ZWIND(:,:)>0.) +#ifndef MNH_BITREP PSFU(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZUA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB) PSFV(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZVA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB) +#else + PSFU(:,:) = - SQRT(BR_P2(ZSFU)+BR_P2(ZSFV)) * ZUA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB) + PSFV(:,:) = - SQRT(BR_P2(ZSFU)+BR_P2(ZSFV)) * ZVA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB) +#endif END WHERE ! !* conversion from H (W/m2) to w'Theta' diff --git a/src/MNH/mode_prandtl.f90 b/src/MNH/mode_prandtl.f90 index d1665a314a0f5be73bf08f267c3a5e365b6f913d..32b9fb4ef03e2643efcca2fed9197c7e6f4a550a 100644 --- a/src/MNH/mode_prandtl.f90 +++ b/src/MNH/mode_prandtl.f90 @@ -17,6 +17,9 @@ ! USE MODD_CTURB, ONLY : XCTV, XCSHF, XCTD, XPHI_LIM, XCPR3, XCPR4, XCPR5 USE MODD_PARAMETERS, ONLY : JPVEXT_TURB +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! USE MODI_SHUMAN IMPLICIT NONE @@ -80,7 +83,11 @@ IF (HTURBDIM=='3DIM') THEN !* 3DIM case IF (OUSERV) THEN ZW1(:,:,:) = 1. + 1.5* (PREDTH1(:,:,:)+PREDR1(:,:,:)) + & +#ifndef MNH_BITREP ( 0.5 * (PREDTH1(:,:,:)**2+PREDR1(:,:,:)**2) & +#else + ( 0.5 * (BR_P2(PREDTH1(:,:,:))+BR_P2(PREDR1(:,:,:))) & +#endif + PREDTH1(:,:,:) * PREDR1(:,:,:) & ) ZW2(:,:,:) = 0.5 * (PRED2TH3(:,:,:)-PRED2R3(:,:,:)) @@ -91,7 +98,11 @@ IF (HTURBDIM=='3DIM') THEN ) / ZW1(:,:,:) ELSE ZW1(:,:,:) = 1. + 1.5* PREDTH1(:,:,:) + & +#ifndef MNH_BITREP 0.5* PREDTH1(:,:,:)**2 +#else + 0.5* BR_P2(PREDTH1(:,:,:)) +#endif ZW2(:,:,:) = 0.5* PRED2TH3(:,:,:) PPHI3(:,:,:)= 1. - & (PRED2TH3(:,:,:) / PREDTH1(:,:,:) + ZW2(:,:,:)) / ZW1(:,:,:) @@ -588,13 +599,21 @@ IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC PD_M3_WTH_W2TH_O_DDTDZ(:,:,:) = & +#ifndef MNH_BITREP - XCSHF*PKEFF*1.5/MZM(KKA,KKU,KKL,PTKE)/(1.+PREDTH1)**2*XCTV*PBLL_O_E*PETHETA & +#else + - XCSHF*PKEFF*1.5/MZM(KKA,KKU,KKL,PTKE)/BR_P2(1.+PREDTH1)*XCTV*PBLL_O_E*PETHETA & +#endif * (1. - 0.5*PREDR1*(1.+PREDR1)/PD*( 1.+(1.+PREDTH1)*(1.5+PREDR1+PREDTH1)/PD) ) #else CALL MZM_DEVICE(PTKE,ZTMP1_DEVICE) !$acc kernels PD_M3_WTH_W2TH_O_DDTDZ(:,:,:) = & +#ifndef MNH_BITREP - XCSHF*PKEFF*1.5/ZTMP1_DEVICE/(1.+PREDTH1)**2*XCTV*PBLL_O_E*PETHETA & +#else + - XCSHF*PKEFF*1.5/ZTMP1_DEVICE/BR_P2(1.+PREDTH1)*XCTV*PBLL_O_E*PETHETA & +#endif * (1. - 0.5*PREDR1*(1.+PREDR1)/PD*( 1.+(1.+PREDTH1)*(1.5+PREDR1+PREDTH1)/PD) ) #endif ! @@ -751,7 +770,11 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PM3_WTH_WR2(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*PEMOIST**2 & +#else +PM3_WTH_WR2(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*BR_P2(PEMOIST) & +#endif *MZM(KKA,KKU,KKL,PBETA*PLEPS/(PSQRT_TKE*PTKE))/XCTD*PDTDZ/PD #else !$acc kernels @@ -759,7 +782,11 @@ ZTMP2_DEVICE = PBETA*PLEPS/(PSQRT_TKE*PTKE) !$acc end kernels CALL MZM_DEVICE(ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PM3_WTH_WR2(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*PEMOIST**2 & +#else +PM3_WTH_WR2(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*BR_P2(PEMOIST) & +#endif *ZTMP1_DEVICE/XCTD*PDTDZ/PD #endif ! @@ -810,7 +837,11 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PD_M3_WTH_WR2_O_DDTDZ(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*PEMOIST**2 & +#else +PD_M3_WTH_WR2_O_DDTDZ(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*BR_P2(PEMOIST) & +#endif *MZM(KKA,KKU,KKL,PBETA*PLEPS/(PSQRT_TKE*PTKE))/XCTD/PD & * (1. - PREDTH1*(1.5+PREDTH1+PREDR1)/PD) #else @@ -819,7 +850,11 @@ ZTMP2_DEVICE = PBETA*PLEPS/(PSQRT_TKE*PTKE) !$acc end kernels CALL MZM_DEVICE(ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PD_M3_WTH_WR2_O_DDTDZ(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*PEMOIST**2 & +#else +PD_M3_WTH_WR2_O_DDTDZ(:,:,:) = - XCSHF*PKEFF*0.25*PBLL_O_E*XCTV*BR_P2(PEMOIST) & +#endif *ZTMP1_DEVICE/XCTD/PD & * (1. - PREDTH1*(1.5+PREDTH1+PREDR1)/PD) #endif @@ -1021,16 +1056,28 @@ IKE = SIZE(PD,3)-JPVEXT_TURB IF (OUSERV) THEN PD_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM*PLEPS/PTKE*XCTV * MZF(KKA,KKU,KKL, & (1.-0.5*PREDR1*(1.+PREDR1)/PD)*(1.-(1.5+PREDTH1+PREDR1)* & +#ifndef MNH_BITREP PREDTH1*(1.+PREDTH1)/PD ) / (1.+PREDTH1)**2 ) +#else + PREDTH1*(1.+PREDTH1)/PD ) / BR_P2(1.+PREDTH1) ) +#endif ELSE +#ifndef MNH_BITREP PD_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM*PLEPS/PTKE*XCTV * MZF(KKA,KKU,KKL,1./(1.+PREDTH1)**2) +#else + PD_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM*PLEPS/PTKE*XCTV * MZF(KKA,KKU,KKL,1./BR_P2(1.+PREDTH1)) +#endif END IF #else IF (OUSERV) THEN !$acc kernels ZTMP2_DEVICE = (1.-0.5*PREDR1*(1.+PREDR1)/PD)*(1.-(1.5+PREDTH1+PREDR1)* & +#ifndef MNH_BITREP PREDTH1*(1.+PREDTH1)/PD ) / (1.+PREDTH1)**2 +#else + PREDTH1*(1.+PREDTH1)/PD ) / BR_P2(1.+PREDTH1) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE, ZTMP1_DEVICE) !$acc kernels @@ -1039,7 +1086,11 @@ CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE, ZTMP1_DEVICE) ELSE !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = 1./(1.+PREDTH1)**2 +#else +ZTMP2_DEVICE = 1./BR_P2(1.+PREDTH1) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels @@ -1092,10 +1143,18 @@ IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC PM3_TH2_WTH2(:,:,:) = PLEPS*0.5/XCTD/PSQRT_TKE & +#ifndef MNH_BITREP * MZF(KKA,KKU,KKL, (1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD ) +#else + * MZF(KKA,KKU,KKL, (1.+0.5*PREDTH1+1.5*PREDR1+0.5*BR_P2(PREDR1))/PD ) +#endif #else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = (1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD +#else +ZTMP2_DEVICE = (1.+0.5*PREDTH1+1.5*PREDR1+0.5*BR_P2(PREDR1))/PD +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL, ZTMP2_DEVICE ,ZTMP1_DEVICE) !$acc kernels @@ -1148,12 +1207,20 @@ IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC PD_M3_TH2_WTH2_O_DDTDZ(:,:,:) = PLEPS*0.5/XCTD/PSQRT_TKE*XCTV & * MZF(KKA,KKU,KKL, PBLL_O_E*PETHETA* (0.5/PD & +#ifndef MNH_BITREP - (1.5+PREDTH1+PREDR1)*(1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD**2 & +#else + - (1.5+PREDTH1+PREDR1)*(1.+0.5*PREDTH1+1.5*PREDR1+0.5*BR_P2(PREDR1))/BR_P2(PD) & +#endif ) ) #else !$acc kernels ZTMP2_DEVICE = PBLL_O_E*PETHETA* (0.5/PD & +#ifndef MNH_BITREP - (1.5+PREDTH1+PREDR1)*(1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD**2 & +#else + - (1.5+PREDTH1+PREDR1)*(1.+0.5*PREDTH1+1.5*PREDR1+0.5*BR_P2(PREDR1))/BR_P2(PD) & +#endif ) !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL, ZTMP2_DEVICE ,ZTMP1_DEVICE) @@ -1205,14 +1272,26 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PM3_TH2_W2R(:,:,:) = 0.75*XCTV**2*MZF(KKA,KKU,KKL,PBLL_O_E*PEMOIST/PD*PDTDZ**2)*PLM*PLEPS/PTKE #else +PM3_TH2_W2R(:,:,:) = 0.75*BR_P2(XCTV)*MZF(KKA,KKU,KKL,PBLL_O_E*PEMOIST/PD*BR_P2(PDTDZ))*PLM*PLEPS/PTKE +#endif +#else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = PBLL_O_E*PEMOIST/PD*PDTDZ**2 +#else +ZTMP2_DEVICE = PBLL_O_E*PEMOIST/PD*BR_P2(PDTDZ) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PM3_TH2_W2R(:,:,:) = 0.75*XCTV**2*ZTMP1_DEVICE*PLM*PLEPS/PTKE +#else +PM3_TH2_W2R(:,:,:) = 0.75*BR_P2(XCTV)*ZTMP1_DEVICE*PLM*PLEPS/PTKE +#endif #endif ! PM3_TH2_W2R(:,:,IKB-1)=PM3_TH2_W2R(:,:,IKB) @@ -1261,7 +1340,11 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PD_M3_TH2_W2R_O_DDTDZ(:,:,:) = 0.75*XCTV**2*PLM*PLEPS/PTKE & +#else +PD_M3_TH2_W2R_O_DDTDZ(:,:,:) = 0.75*BR_P2(XCTV)*PLM*PLEPS/PTKE & +#endif * MZF(KKA,KKU,KKL, PBLL_O_E*PEMOIST/PD*PDTDZ*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) ) #else !$acc kernels @@ -1269,7 +1352,11 @@ ZTMP2_DEVICE = PBLL_O_E*PEMOIST/PD*PDTDZ*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PD_M3_TH2_W2R_O_DDTDZ(:,:,:) = 0.75*XCTV**2*PLM*PLEPS/PTKE & +#else +PD_M3_TH2_W2R_O_DDTDZ(:,:,:) = 0.75*BR_P2(XCTV)*PLM*PLEPS/PTKE & +#endif * ZTMP1_DEVICE #endif ! @@ -1316,14 +1403,26 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PM3_TH2_WR2(:,:,:) = 0.25*XCTV**2*MZF(KKA,KKU,KKL,(PBLL_O_E*PEMOIST*PDTDZ)**2/PD)*PLEPS/PSQRT_TKE/XCTD #else +PM3_TH2_WR2(:,:,:) = 0.25*BR_P2(XCTV)*MZF(KKA,KKU,KKL,BR_P2(PBLL_O_E*PEMOIST*PDTDZ)/PD)*PLEPS/PSQRT_TKE/XCTD +#endif +#else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = (PBLL_O_E*PEMOIST*PDTDZ)**2/PD +#else +ZTMP2_DEVICE = BR_P2(PBLL_O_E*PEMOIST*PDTDZ)/PD +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PM3_TH2_WR2(:,:,:) = 0.25*XCTV**2*ZTMP1_DEVICE*PLEPS/PSQRT_TKE/XCTD +#else +PM3_TH2_WR2(:,:,:) = 0.25*BR_P2(XCTV)*ZTMP1_DEVICE*PLEPS/PSQRT_TKE/XCTD +#endif #endif ! PM3_TH2_WR2(:,:,IKB-1)=PM3_TH2_WR2(:,:,IKB) @@ -1371,15 +1470,28 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PD_M3_TH2_WR2_O_DDTDZ(:,:,:) = 0.25*XCTV**2*PLEPS/PSQRT_TKE/XCTD & * MZF(KKA,KKU,KKL, (PBLL_O_E*PEMOIST)**2*PDTDZ/PD*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) ) #else +PD_M3_TH2_WR2_O_DDTDZ(:,:,:) = 0.25*BR_P2(XCTV)*PLEPS/PSQRT_TKE/XCTD & + * MZF(KKA,KKU,KKL, BR_P2(PBLL_O_E*PEMOIST)*PDTDZ/PD*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) ) +#endif +#else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = (PBLL_O_E*PEMOIST)**2*PDTDZ/PD*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) +#else +ZTMP2_DEVICE = BR_P2(PBLL_O_E*PEMOIST)*PDTDZ/PD*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PD_M3_TH2_WR2_O_DDTDZ(:,:,:) = 0.25*XCTV**2*PLEPS/PSQRT_TKE/XCTD & +#else +PD_M3_TH2_WR2_O_DDTDZ(:,:,:) = 0.25*BR_P2(XCTV)*PLEPS/PSQRT_TKE/XCTD & +#endif * ZTMP1_DEVICE #endif ! @@ -1703,15 +1815,28 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PD_M3_THR_WTH2_O_DDTDZ(:,:,:) = - 0.25*PLEPS/PSQRT_TKE/XCTD*XCTV**2 & * MZF(KKA,KKU,KKL, -(1.+PREDR1)*(PBLL_O_E*PETHETA/PD)**2*PDRDZ*(1.5+PREDTH1+PREDR1) ) #else +PD_M3_THR_WTH2_O_DDTDZ(:,:,:) = - 0.25*PLEPS/PSQRT_TKE/XCTD*BR_P2(XCTV) & + * MZF(KKA,KKU,KKL, -(1.+PREDR1)*BR_P2(PBLL_O_E*PETHETA/PD)*PDRDZ*(1.5+PREDTH1+PREDR1) ) +#endif +#else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = -(1.+PREDR1)*(PBLL_O_E*PETHETA/PD)**2*PDRDZ*(1.5+PREDTH1+PREDR1) +#else +ZTMP2_DEVICE = -(1.+PREDR1)*BR_P2(PBLL_O_E*PETHETA/PD)*PDRDZ*(1.5+PREDTH1+PREDR1) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PD_M3_THR_WTH2_O_DDTDZ(:,:,:) = - 0.25*PLEPS/PSQRT_TKE/XCTD*XCTV**2 * ZTMP1_DEVICE +#else +PD_M3_THR_WTH2_O_DDTDZ(:,:,:) = - 0.25*PLEPS/PSQRT_TKE/XCTD*BR_P2(XCTV) * ZTMP1_DEVICE +#endif #endif ! PD_M3_THR_WTH2_O_DDTDZ(:,:,IKB-1)=PD_M3_THR_WTH2_O_DDTDZ(:,:,IKB) @@ -1872,15 +1997,28 @@ IKB = 1+JPVEXT_TURB IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC +#ifndef MNH_BITREP PD_M3_THR_W2TH_O_DDTDZ(:,:,:) = - 0.75*PLM*PLEPS/PTKE * XCTV**2 & * MZF(KKA,KKU,KKL, -PETHETA*PBLL_O_E*(1.+PREDR1)*PDRDZ*(1.5+PREDTH1+PREDR1)/PD**2 ) #else +PD_M3_THR_W2TH_O_DDTDZ(:,:,:) = - 0.75*PLM*PLEPS/PTKE * BR_P2(XCTV) & + * MZF(KKA,KKU,KKL, -PETHETA*PBLL_O_E*(1.+PREDR1)*PDRDZ*(1.5+PREDTH1+PREDR1)/BR_P2(PD) ) +#endif +#else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = -PETHETA*PBLL_O_E*(1.+PREDR1)*PDRDZ*(1.5+PREDTH1+PREDR1)/PD**2 +#else +ZTMP2_DEVICE = -PETHETA*PBLL_O_E*(1.+PREDR1)*PDRDZ*(1.5+PREDTH1+PREDR1)/BR_P2(PD) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP PD_M3_THR_W2TH_O_DDTDZ(:,:,:) = - 0.75*PLM*PLEPS/PTKE * XCTV**2 * ZTMP1_DEVICE +#else +PD_M3_THR_W2TH_O_DDTDZ(:,:,:) = - 0.75*PLM*PLEPS/PTKE * BR_P2(XCTV) * ZTMP1_DEVICE +#endif #endif ! PD_M3_THR_W2TH_O_DDTDZ(:,:,IKB-1)=PD_M3_THR_W2TH_O_DDTDZ(:,:,IKB) @@ -1927,12 +2065,20 @@ IKE = SIZE(PD,3)-JPVEXT_TURB #ifndef _OPENACC PD_M3_THR_W2TH_O_DDRDZ(:,:,:) = - 0.75*PLM*PLEPS/PTKE * XCTV & +#ifndef MNH_BITREP * MZF(KKA,KKU,KKL, -(1.+PREDR1)*PREDR1*(1.5+PREDTH1+PREDR1)/PD**2 & +#else + * MZF(KKA,KKU,KKL, -(1.+PREDR1)*PREDR1*(1.5+PREDTH1+PREDR1)/BR_P2(PD) & +#endif +(1.+2.*PREDR1)/PD & ) #else !$acc kernels +#ifndef MNH_BITREP ZTMP2_DEVICE = -(1.+PREDR1)*PREDR1*(1.5+PREDTH1+PREDR1)/PD**2 & +#else +ZTMP2_DEVICE = -(1.+PREDR1)*PREDR1*(1.5+PREDTH1+PREDR1)/BR_P2(PD) & +#endif +(1.+2.*PREDR1)/PD !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP2_DEVICE,ZTMP1_DEVICE) diff --git a/src/MNH/ppm.f90 b/src/MNH/ppm.f90 index 591a141a8226360e1fa28792b486298f5841b163..917bca26486d8d3b65596da8649d2d3bf7e511de 100644 --- a/src/MNH/ppm.f90 +++ b/src/MNH/ppm.f90 @@ -443,6 +443,9 @@ USE MODE_MPPDB #ifdef _OPENACC USE MODD_PARAMETERS, ONLY : JPHEXT #endif +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -631,11 +634,19 @@ PRINT *,'OPENACC: ppm::PPM_01_X CYCL/WALL boundaries not yet implemented' ZQL = PSRC ZQR = PSRC ZQ6 = 0.0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ < -BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQL0 - PSRC) ZQR = ZQL0 - ZQ6 ZQL = ZQL0 - ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#ifndef MNH_BITREP + ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ > BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQR0 - PSRC) ZQL = ZQR0 - ZQ6 ZQR = ZQR0 @@ -799,11 +810,19 @@ CASE('OPEN') ZQL = PSRC ZQR = PSRC ZQ6 = 0.0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ < -BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQL0 - PSRC) ZQR = ZQL0 - ZQ6 ZQL = ZQL0 - ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#ifndef MNH_BITREP + ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ > BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQR0 - PSRC) ZQL = ZQR0 - ZQ6 ZQR = ZQR0 @@ -835,11 +854,19 @@ DO K=1,IKU ZQL(I,J,K) = PSRC(I,J,K) ZQR(I,J,K) = PSRC(I,J,K) ZQ6(I,J,K) = 0.0 +#ifndef MNH_BITREP ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -(ZDQ(I,J,K))**2 ) THEN +#else + ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K)) ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K) ZQL(I,J,K) = ZQL0(I,J,K) +#ifndef MNH_BITREP ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > (ZDQ(I,J,K))**2 ) THEN +#else + ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K)) ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K) ZQR(I,J,K) = ZQR0(I,J,K) @@ -1006,6 +1033,9 @@ END FUNCTION PPM_01_X USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU +#ifdef MNH_BITREP +USE MODI_BITREP +#endif IMPLICIT NONE ! @@ -1277,11 +1307,19 @@ PRINT *,'OPENACC: ppm::PPM_01_Y CYCL/WALL boundaries not yet implemented' ZQL = PSRC ZQR = PSRC ZQ6 = 0.0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ < -BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQL0 - PSRC) ZQR = ZQL0 - ZQ6 ZQL = ZQL0 - ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#ifndef MNH_BITREP + ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ > BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQR0 - PSRC) ZQL = ZQR0 - ZQ6 ZQR = ZQR0 @@ -1422,11 +1460,19 @@ CALL GET_HALO_D(ZQL0,HDIR="01_Y") ZQL = PSRC ZQR = PSRC ZQ6 = 0.0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ < -BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQL0 - PSRC) ZQR = ZQL0 - ZQ6 ZQL = ZQL0 - ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#ifndef MNH_BITREP + ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) +#else + ELSEWHERE ( ZQ60*ZDQ > BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQR0 - PSRC) ZQL = ZQR0 - ZQ6 ZQR = ZQR0 @@ -1462,12 +1508,20 @@ CALL GET_HALO_D(ZQL0,HDIR="01_Y") ZQR(I,J,K) = PSRC(I,J,K) ZQ6(I,J,K) = 0.0 END IF +#ifndef MNH_BITREP IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -(ZDQ(I,J,K))**2 ) ) THEN +#else + IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K)) ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K) ZQL(I,J,K) = ZQL0(I,J,K) END IF +#ifndef MNH_BITREP IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > (ZDQ(I,J,K))**2 ) ) THEN +#else + IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K)) ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K) ZQR(I,J,K) = ZQR0(I,J,K) @@ -1650,7 +1704,9 @@ END FUNCTION PPM_01_Y USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU - +#ifdef MNH_BITREP +USE MODI_BITREP +#endif IMPLICIT NONE ! @@ -1852,12 +1908,19 @@ WHERE ( ZDMQ == 0.0 ) ZQL = PSRC ZQR = PSRC ZQ6 = 0.0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ < -(ZDQ)**2 ) +#else +ELSEWHERE ( ZQ60*ZDQ < -BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQL0 - PSRC) ZQR = ZQL0 - ZQ6 ZQL = ZQL0 +#ifndef MNH_BITREP ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) #else +ELSEWHERE ( ZQ60*ZDQ > BR_P2(ZDQ) ) +#endif ZQ6 = 3.0*(ZQR0 - PSRC) ZQL = ZQR0 - ZQ6 ZQR = ZQR0 @@ -1884,13 +1947,21 @@ ZDQ = ZQR - ZQL !!$ ZQR = PSRC !!$ ZQ6 = 0.0 !!$END WHERE +!!#ifndef MNH_BITREP !!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -ZDQ**2 ) ) +!!#else +!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -BR_P2(ZDQ) ) ) +!!#endif !!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -BR_P2(ZDQ) ) ) !!$ ZQ6 = 3.0*(ZQL0 - PSRC) !!$ ZQR = ZQL0 - ZQ6 !!$ ZQL = ZQL0 !!$END WHERE +!!#ifndef MNH_BITREP !!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > ZDQ**2 ) ) +!!#else +!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > BR_P2(ZDQ) ) ) +!!#endif !!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > BR_P2(ZDQ) ) ) !!$ ZQ6 = 3.0*(ZQR0 - PSRC) !!$ ZQL = ZQR0 - ZQ6 @@ -1925,12 +1996,20 @@ ZDQ = ZQR - ZQL ZQR(I,J,K) = PSRC(I,J,K) ZQ6(I,J,K) = 0.0 END IF +#ifndef MNH_BITREP IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)**2 ) ) THEN +#else + IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K)) ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K) ZQL(I,J,K) = ZQL0(I,J,K) END IF +#ifndef MNH_BITREP IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)**2 ) ) THEN +#else + IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) ) THEN +#endif ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K)) ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K) ZQR(I,J,K) = ZQR0(I,J,K) diff --git a/src/MNH/prandtl.f90 b/src/MNH/prandtl.f90 index 4eec1343f428ccd855a74b1fa0a3df9b80e3cd01..7349a546e9662a62b8482fc63385ac9b8579cb87 100644 --- a/src/MNH/prandtl.f90 +++ b/src/MNH/prandtl.f90 @@ -218,6 +218,9 @@ USE MODI_SHUMAN USE MODI_SHUMAN_DEVICE #endif USE MODE_FMWRIT +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -454,11 +457,19 @@ PRINT *,'OPENACC: prandtl::1DIM not yet tested' #endif ! !$acc kernels async +#ifndef MNH_BITREP PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +#else + PRED2TH3(:,:,:) = BR_P2(PREDTH1(:,:,:)) +#endif !$acc end kernels ! !$acc kernels async +#ifndef MNH_BITREP PRED2R3(:,:,:) = PREDR1(:,:,:) **2 +#else + PRED2R3(:,:,:) = BR_P2(PREDR1(:,:,:)) +#endif !$acc end kernels ! !$acc kernels async @@ -473,39 +484,69 @@ ELSE IF (L2D) THEN ! 3D case in a 2D model PRINT *,'OPENACC: prandtl::L2D=.T. and KRR/=0 not yet tested' #endif #ifndef _OPENACC +#ifndef MNH_BITREP PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2+(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)**2 ) +#else + PRED2TH3(:,:,:)= BR_P2(PREDTH1(:,:,:))+BR_P2(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) ) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)) ) +#endif PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE ) !$acc kernels async +#ifndef MNH_BITREP PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2+(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * ZTMP2_DEVICE +#else + PRED2TH3(:,:,:)= BR_P2(PREDTH1(:,:,:))+BR_P2(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) ) * ZTMP2_DEVICE +#endif PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) !$acc end kernels #endif ! #ifndef _OPENACC +#ifndef MNH_BITREP PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 ) +#else + PRED2R3(:,:,:)= BR_P2(PREDR1(:,:,:)) + BR_P2(XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:)) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)) ) +#endif PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL) #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE ) !$acc kernels async +#ifndef MNH_BITREP PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * ZTMP2_DEVICE +#else + PRED2R3(:,:,:)= BR_P2(PREDR1(:,:,:)) + BR_P2(XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:)) * ZTMP2_DEVICE +#endif PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL) !$acc end kernels #endif ! #ifndef _OPENACC +#ifndef MNH_BITREP PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + XCTV**2*PBLL_O_E(:,:,:)**2 * & +#else + PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & +#endif PEMOIST(:,:,:) * PETHETA(:,:,:) * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)* & GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)) @@ -518,7 +559,11 @@ PRINT *,'OPENACC: prandtl::L2D=.T. and KRR/=0 not yet tested' !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels async +#ifndef MNH_BITREP PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + XCTV**2*PBLL_O_E(:,:,:)**2 * & +#else + PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & +#endif PEMOIST(:,:,:) * PETHETA(:,:,:) * ZTMP2_DEVICE PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+KKL) !$acc end kernels @@ -530,17 +575,30 @@ PRINT *,'OPENACC: prandtl::L2D=.T. and KRR/=0 not yet tested' PRINT *,'OPENACC: prandtl::L2D=.T. and KRR=0 not yet tested' #endif #ifndef _OPENACC +#ifndef MNH_BITREP PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 + XCTV**2*PBLL_O_E(:,:,:)**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)**2 ) +#else + PRED2TH3(:,:,:) = BR_P2'PREDTH1(:,:,:)) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)) ) +#endif PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 + XCTV**2*PBLL_O_E(:,:,:)**2 * ZTMP2_DEVICE +#else + PRED2TH3(:,:,:) = BR_P2(PREDTH1(:,:,:)) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * ZTMP2_DEVICE +#endif !$acc end kernels !$acc kernels async PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) @@ -562,37 +620,65 @@ ELSE ! 3D case in a 3D model ! IF (KRR /= 0) THEN ! moist 3D case #ifndef _OPENACC +#ifndef MNH_BITREP PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2 + ( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)**2 & + GY_M_M(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY)**2 ) +#else + PRED2TH3(:,:,:)= BR_P2(PREDTH1(:,:,:)) + BR_P2( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) ) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)) & + + BR_P2(GY_M_M(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY)) ) +#endif PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) CALL GY_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 + ZTMP2_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) + BR_P2(ZTMP2_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2 + ( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * ZTMP2_DEVICE +#else + PRED2TH3(:,:,:)= BR_P2(PREDTH1(:,:,:)) + BR_P2( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) ) * ZTMP2_DEVICE +#endif PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL) !$acc end kernels #endif ! #ifndef _OPENACC +#ifndef MNH_BITREP PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 + & GY_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDYY,PDZZ,PDZY)**2 ) +#else + PRED2R3(:,:,:)= BR_P2(PREDR1(:,:,:)) + BR_P2(XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:)) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)) + & + BR_P2(GY_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDYY,PDZZ,PDZY)) ) +#endif PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL) #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE) CALL GY_M_M_DEVICE(KKA,KKU,KKL,PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 + ZTMP2_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) + BR_P2(ZTMP2_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * ZTMP2_DEVICE +#else + PRED2R3(:,:,:)= BR_P2(PREDR1(:,:,:)) + BR_P2(XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:)) * ZTMP2_DEVICE +#endif !$acc end kernels !$acc kernels async PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL) @@ -600,7 +686,11 @@ ELSE ! 3D case in a 3D model #endif ! #ifndef _OPENACC +#ifndef MNH_BITREP PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + XCTV**2*PBLL_O_E(:,:,:)**2 * & +#else + PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & +#endif PEMOIST(:,:,:) * PETHETA(:,:,:) * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PRM(:,:,:,1),PDXX,PDZZ,PDZX)* & GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)+ & @@ -617,7 +707,11 @@ ELSE ! 3D case in a 3D model !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + XCTV**2*PBLL_O_E(:,:,:)**2 * & +#else + PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & +#endif PEMOIST(:,:,:) * PETHETA(:,:,:) * ZTMP2_DEVICE !$acc end kernels !$acc kernels async @@ -631,18 +725,32 @@ ELSE ! 3D case in a 3D model PRINT *,'OPENACC: prandtl::L2D=.F. and KRR=0 not yet tested' #endif #ifndef _OPENACC +#ifndef MNH_BITREP PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 + XCTV**2*PBLL_O_E(:,:,:)**2 * & MZM(KKA,KKU,KKL, GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)**2 & + GY_M_M(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY)**2 ) +#else + PRED2TH3(:,:,:) = BR_P2(PREDTH1(:,:,:)) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * & + MZM(KKA,KKU,KKL, BR_P2(GX_M_M(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX)) & + + BR_P2(GY_M_M(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY)) ) +#endif #else CALL GX_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) CALL GY_M_M_DEVICE(KKA,KKU,KKL,PTHLM,PDYY,PDZZ,PDZY,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = ZTMP1_DEVICE**2 + ZTMP2_DEVICE**2 +#else + ZTMP1_DEVICE = BR_P2(ZTMP1_DEVICE) + BR_P2(ZTMP2_DEVICE) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels +#ifndef MNH_BITREP PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 + XCTV**2*PBLL_O_E(:,:,:)**2 * ZTMP2_DEVICE +#else + PRED2TH3(:,:,:) = BR_P2(PREDTH1(:,:,:)) + BR_P2(XCTV)*BR_P2(PBLL_O_E(:,:,:)) * ZTMP2_DEVICE +#endif !$acc end kernels #endif !$acc kernels async @@ -685,15 +793,27 @@ ELSE IF (L2D) THEN ! 3D case in a 2D model ! #ifndef _OPENACC DO JSV=1,ISV +#ifndef MNH_BITREP IF (KRR /= 0) THEN ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA ELSE ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2) END IF +#else + IF (KRR /= 0) THEN + ZW1 = MZM(KKA,KKU,KKL, BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM) ) *PETHETA + ELSE + ZW1 = MZM(KKA,KKU,KKL, BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM)) + END IF +#endif #else DO JSV=1,ISV !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 +#else + ZTMP1_DEVICE = BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZW1) IF (KRR /= 0) THEN @@ -751,15 +871,27 @@ ELSE ! 3D case in a 3D model ! #ifndef _OPENACC DO JSV=1,ISV +#ifndef MNH_BITREP IF (KRR /= 0) THEN ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA ELSE ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2) END IF +#else + IF (KRR /= 0) THEN + ZW1 = MZM(KKA,KKU,KKL, BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM) ) *PETHETA + ELSE + ZW1 = MZM(KKA,KKU,KKL, BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM)) + END IF +#endif #else DO JSV=1,ISV !$acc kernels +#ifndef MNH_BITREP ZTMP1_DEVICE = (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 +#else + ZTMP1_DEVICE = BR_P2(XG / PTHVREF * PLM * PLEPS / PTKEM) +#endif !$acc end kernels CALL MZM_DEVICE(ZTMP1_DEVICE,ZW1) IF (KRR /= 0) THEN diff --git a/src/MNH/rain_ice.f90 b/src/MNH/rain_ice.f90 index ace041532958db20731d1084be8cb780f3c862bf..891a9856dc2d4c8a5b381aa3af323b54608c8752 100644 --- a/src/MNH/rain_ice.f90 +++ b/src/MNH/rain_ice.f90 @@ -2298,9 +2298,9 @@ IMPLICIT NONE ZCJ(:) = XSCFAC * ZRHODREF(:)**0.3 / SQRT( 1.718E-5+0.0049E-5*(ZZT(:)-XTT) ) #else ZAI(:) = BR_EXP( XALPI - XBETAI/ZZT(:) - XGAMI*BR_LOG(ZZT(:) ) ) ! es_i - ZAI(:) = ( XLSTT + (XCPV-XCI)*(ZZT(:)-XTT) )**2 / (ZKA(:)*XRV*ZZT(:)**2) & + ZAI(:) = BR_P2( XLSTT + (XCPV-XCI)*(ZZT(:)-XTT) ) / (ZKA(:)*XRV*BR_P2(ZZT(:))) & + ( XRV*ZZT(:) ) / (ZDV(:)*ZAI(:)) - ZCJ(:) = XSCFAC * BR_POW(ZRHODREF(:),0.3) / SQRT( 1.718E-5+0.0049E-5*(ZZT(:)-XTT) ) + ZCJ(:) = XSCFAC * BR_POW(ZRHODREF(:),0.3) / BR_POW( 1.718E-5+0.0049E-5*(ZZT(:)-XTT) , 0.5) #endif ! !* 3.4.2 compute the riming-conversion of r_c for r_i production: RCAUTI @@ -2500,7 +2500,7 @@ REAL :: ZCRIAUTC ! Critical cloud mixing ratio ZZW(JL) = MIN( ZRCS(JL) , XTIMAUTC*( ZRCT(JL)+ZSIGMA_RC(JL)-ZCRIAUTC )**2 & /( 4. * ZSIGMA_RC(JL) ) ) #else - ZZW(JL) = MIN( ZRCS(JL) , XTIMAUTC*( ZRCT(JL)+ZSIGMA_RC(JL)-ZCRIAUTC )**2 & + ZZW(JL) = MIN( ZRCS(JL) , XTIMAUTC*BR_P2( ZRCT(JL)+ZSIGMA_RC(JL)-ZCRIAUTC ) & /( 4. * ZSIGMA_RC(JL) ) ) #endif ENDIF @@ -2574,7 +2574,7 @@ REAL :: ZCRIAUTC ! Critical cloud mixing ratio ZZW(:) = MIN( ZRRS(:),( MAX( 0.0,ZUSW(:) )/(ZRHODREF(:)*ZZW(:)) ) * & ( X0EVAR*ZLBDAR(:)**XEX0EVAR+X1EVAR*ZCJ(:)*ZLBDAR(:)**XEX1EVAR ) ) #else - ZZW(:) = ( XLVTT+(XCPV-XCL)*(ZZT(:)-XTT) )**2 / ( ZKA(:)*XRV*ZZT(:)**2 ) & + ZZW(:) = BR_P2( XLVTT+(XCPV-XCL)*(ZZT(:)-XTT) ) / ( ZKA(:)*XRV*BR_P2(ZZT(:)) ) & + ( XRV*ZZT(:) ) / ( ZDV(:)*ZZW(:) ) ZZW(:) = MIN( ZRRS(:),( MAX( 0.0,ZUSW(:) )/(ZRHODREF(:)*ZZW(:)) ) * & ( X0EVAR*BR_POW(ZLBDAR(:),XEX0EVAR)+X1EVAR*ZCJ(:)*BR_POW(ZLBDAR(:),XEX1EVAR) ) ) @@ -2883,9 +2883,9 @@ REAL,DIMENSION(SIZE(ZZW1,1)) :: ZTMP #else ZZW1(:,2) = & !! coef of RRACCS XFRACCSS*( BR_POW(ZLBDAS(:),XCXS) )*( BR_POW(ZRHODREF(:),-XCEXVT-1.) ) & - *( XLBRACCS1/((ZLBDAS(:)**2) ) + & + *( XLBRACCS1/((BR_P2(ZLBDAS(:))) ) + & XLBRACCS2/( ZLBDAS(:) * ZLBDAR(:) ) + & - XLBRACCS3/( (ZLBDAR(:)**2)) )/ZLBDAR(:)**4 + XLBRACCS3/( (BR_P2(ZLBDAR(:)))) )/BR_POW(ZLBDAR(:),4.0) #endif ZZW1(:,4) = MIN( ZRRS(:),ZZW1(:,2)*ZZW(:) ) ! RRACCSS ZRRS(:) = ZRRS(:) - ZZW1(:,4) @@ -2964,9 +2964,9 @@ REAL,DIMENSION(SIZE(ZZW1,1)) :: ZTMP #else ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)* & ! RSACCRG ( BR_POW(ZLBDAS(:),XCXS-XBS) )*( BR_POW(ZRHODREF(:),-XCEXVT-1.) ) & - *( XLBSACCR1/((ZLBDAR(:)**2) ) + & + *( XLBSACCR1/((BR_P2(ZLBDAR(:))) ) + & XLBSACCR2/( ZLBDAR(:) * ZLBDAS(:) ) + & - XLBSACCR3/( (ZLBDAS(:)**2)) )/ZLBDAR(:) ) + XLBSACCR3/( (BR_P2(ZLBDAS(:)))) )/ZLBDAR(:) ) #endif ZRRS(:) = ZRRS(:) - ZZW1(:,2) ZRSS(:) = ZRSS(:) - ZZW1(:,3) @@ -3235,9 +3235,9 @@ INTEGER,DIMENSION(:),ALLOCATABLE :: I1 * BR_EXP( XCOLEXSG*(ZZT(:)-XTT) ) & *( BR_POW(ZLBDAS(:),XCXS-XBS) )*( BR_POW(ZLBDAG(:),XCXG) ) & *( BR_POW(ZRHODREF(:),-XCEXVT-1.) ) & - *( XLBSDRYG1/( ZLBDAG(:)**2 ) + & + *( XLBSDRYG1/( BR_P2(ZLBDAG(:)) ) + & XLBSDRYG2/( ZLBDAG(:) * ZLBDAS(:) ) + & - XLBSDRYG3/( ZLBDAS(:)**2) ) ) + XLBSDRYG3/( BR_P2(ZLBDAS(:))) ) ) #endif END WHERE !$acc end kernels @@ -3341,9 +3341,9 @@ INTEGER,DIMENSION(:),ALLOCATABLE :: I1 ZZW1(:,4) = MIN( ZRRS(:),XFRDRYG*ZZW(:) & ! RRDRYG *( BR_POW(ZLBDAR(:),-4.) )*( BR_POW(ZLBDAG(:),XCXG) ) & *( BR_POW(ZRHODREF(:),-XCEXVT-1.) ) & - *( XLBRDRYG1/( ZLBDAG(:)**2 ) + & + *( XLBRDRYG1/( BR_P2(ZLBDAG(:)) ) + & XLBRDRYG2/( ZLBDAG(:) * ZLBDAR(:) ) + & - XLBRDRYG3/( ZLBDAR(:)**2) ) ) + XLBRDRYG3/( BR_P2(ZLBDAR(:))) ) ) #endif END WHERE !$acc end kernels @@ -3677,12 +3677,21 @@ CALL ABORT ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 ) ! WHERE( GWET(:) ) +#ifndef MNH_BITREP ZZW1(:,3) = MIN( ZRSS(:),XFSWETH*ZZW(:) & ! RSWETH *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAH(:)**XCXH ) & *( ZRHODREF(:)**(-XCEXVT-1.) ) & *( XLBSWETH1/( ZLBDAH(:)**2 ) + & XLBSWETH2/( ZLBDAH(:) * ZLBDAS(:) ) + & XLBSWETH3/( ZLBDAS(:)**2) ) ) +#else + ZZW1(:,3) = MIN( ZRSS(:),XFSWETH*ZZW(:) & ! RSWETH + *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAH(:)**XCXH ) & + *( ZRHODREF(:)**(-XCEXVT-1.) ) & + *( XLBSWETH1/( BR_P2(ZLBDAH(:)) ) + & + XLBSWETH2/( ZLBDAH(:) * ZLBDAS(:) ) + & + XLBSWETH3/( BR_P2(ZLBDAS(:))) ) ) +#endif END WHERE !$acc end data DEALLOCATE(IVEC2) @@ -3749,12 +3758,21 @@ CALL ABORT ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 ) ! WHERE( GWET(:) ) +#ifndef MNH_BITREP ZZW1(:,5) = MAX(MIN( ZRGS(:),XFGWETH*ZZW(:) & ! RGWETH *( ZLBDAG(:)**(XCXG-XBG) )*( ZLBDAH(:)**XCXH ) & *( ZRHODREF(:)**(-XCEXVT-1.) ) & *( XLBGWETH1/( ZLBDAH(:)**2 ) + & XLBGWETH2/( ZLBDAH(:) * ZLBDAG(:) ) + & XLBGWETH3/( ZLBDAG(:)**2) ) ),0. ) +#else + ZZW1(:,5) = MAX(MIN( ZRGS(:),XFGWETH*ZZW(:) & ! RGWETH + *( ZLBDAG(:)**(XCXG-XBG) )*( ZLBDAH(:)**XCXH ) & + *( ZRHODREF(:)**(-XCEXVT-1.) ) & + *( XLBGWETH1/( BR_P2(ZLBDAH(:)) ) + & + XLBGWETH2/( ZLBDAH(:) * ZLBDAG(:) ) + & + XLBGWETH3/( BR_P2(ZLBDAG(:))) ) ),0. ) +#endif END WHERE !$acc end data DEALLOCATE(IVEC2) diff --git a/src/MNH/rmc01.f90 b/src/MNH/rmc01.f90 index 1b7ea4ed868dc53ce650d20c295b96d626092fba..e49f9d6d38396136251c95729f3bef70c9b5914a 100644 --- a/src/MNH/rmc01.f90 +++ b/src/MNH/rmc01.f90 @@ -80,6 +80,9 @@ USE MODD_CTURB USE MODE_SBL ! USE MODI_SHUMAN +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -235,7 +238,11 @@ END SELECT ! DO JK=1,IKT ZL(:,:,JK) = XKARMAN/SQRT(XALPSBL)/XCMFS & +#ifdef MNH_BITREP * ZZZ(:,:,JK)*PDIRCOSZW(:,:)/(ZPHIM(:,:,JK)**2*SQRT(ZPHIE(:,:,JK))) +#else + * ZZZ(:,:,JK)*PDIRCOSZW(:,:)/(BR_P2(ZPHIM(:,:,JK))*SQRT(ZPHIE(:,:,JK))) +#endif END DO ! PLK(:,:,:)=(1.-ZGAM)*ZL+ZGAM*PLK diff --git a/src/MNH/rotate_wind.f90 b/src/MNH/rotate_wind.f90 index 99de384a0a91cde9d49c4e6f0db5acd1739a7d0a..c71acc211e84a661521ed964dae2542a5786c462 100644 --- a/src/MNH/rotate_wind.f90 +++ b/src/MNH/rotate_wind.f90 @@ -114,6 +114,9 @@ END MODULE MODI_ROTATE_WIND !* 0. DECLARATIONS ! ------------ USE MODD_PARAMETERS +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -232,7 +235,11 @@ DO JJ = IJB,IJE DO JI = IIB,IIE PUSLOPE(JI,JJ) = PCOSSLOPE(JI,JJ) * PDIRCOSZW(JI,JJ) * ZUFIN(JI,JJ) + & PSINSLOPE(JI,JJ) * PDIRCOSZW(JI,JJ) * ZVFIN(JI,JJ) + & +#ifndef MNH_BITREP SQRT(1.-PDIRCOSZW(JI,JJ)**2) * ZWFIN(JI,JJ) +#else + SQRT(1.-BR_P2(PDIRCOSZW(JI,JJ))) * ZWFIN(JI,JJ) +#endif ! PVSLOPE(JI,JJ) =-PSINSLOPE(JI,JJ) * ZUFIN(JI,JJ) + & PCOSSLOPE(JI,JJ) * ZVFIN(JI,JJ) diff --git a/src/MNH/tke_eps_sources.f90 b/src/MNH/tke_eps_sources.f90 index 5d49b954eef2f30620ea35374d89a5ae2ccc2d00..bdcb7413b5f01ec8d4f3b28875ee6f591a1e0466 100644 --- a/src/MNH/tke_eps_sources.f90 +++ b/src/MNH/tke_eps_sources.f90 @@ -207,6 +207,10 @@ USE MODE_ll USE MODD_ARGSLIST_ll, ONLY : LIST_ll ! USE MODI_GET_HALO +#ifdef MNH_BITREP +USE MODI_BITREP +#endif +USE MODE_MPPDB ! IMPLICIT NONE ! @@ -337,6 +341,8 @@ ZSOURCE(:,:,:) = PRTKES(:,:,:) / PRHODJ(:,:,:) + PRTKESM(:,:,:) / PRHODJ(:,:,: - PTKEM(:,:,:) / PTSTEP & + PDP(:,:,:) + PTP(:,:,:) + PTR(:,:,:) - PEXPL * ZFLX(:,:,:) * PTKEM(:,:,:) !$acc end kernels +CALL MPPDB_CHECK3DM("tke_eps_sources::",PRECISION,& + & PRTKES,PRHODJ,PRTKESM,PTKEM,PDP,PTP,PTR,ZFLX) ! !* 2.2 implicit vertical TKE transport ! @@ -346,12 +352,20 @@ ZSOURCE(:,:,:) = PRTKES(:,:,:) / PRHODJ(:,:,:) + PRTKESM(:,:,:) / PRHODJ(:,:,: ! #ifndef _OPENACC ZA(:,:,:) = - PTSTEP * XCET * & +#ifdef MNH_BITREP MZM(KKA,KKU,KKL,ZKEFF) * MZM(KKA,KKU,KKL,PRHODJ) / PDZZ**2 +#else + MZM(KKA,KKU,KKL,ZKEFF) * MZM(KKA,KKU,KKL,PRHODJ) / BR_P2(PDZZ) +#endif #else CALL MZM_DEVICE(ZKEFF, ZTMP1_DEVICE) !Warning: re-used later CALL MZM_DEVICE(PRHODJ,ZTMP2_DEVICE) !Warning: re-used later !$acc kernels +#ifdef MNH_BITREP ZA(:,:,:) = - PTSTEP * XCET * ZTMP1_DEVICE * ZTMP2_DEVICE / PDZZ**2 +#else +ZA(:,:,:) = - PTSTEP * XCET * ZTMP1_DEVICE * ZTMP2_DEVICE / BR_P2(PDZZ) +#endif !$acc end kernels #endif ! diff --git a/src/MNH/tm06.f90 b/src/MNH/tm06.f90 index eb82548d9a8a51fd1c1a803ffdb957db48d26e5a..008d8895f1e2b46e096467f53c1de3b9d40f91a1 100644 --- a/src/MNH/tm06.f90 +++ b/src/MNH/tm06.f90 @@ -76,7 +76,9 @@ END MODULE MODI_TM06 USE MODD_PARAMETERS, ONLY : XUNDEF USE MODD_CST, ONLY : XG USE MODD_PARAMETERS, ONLY : JPVEXT_TURB - +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! ! IMPLICIT NONE @@ -138,6 +140,7 @@ END DO !* w'th'2 ! PMTH2 = 0. +#ifdef MNH_BITREP WHERE(ZZ_O_H < 0.95 .AND. ZZ_O_H/=XUNDEF) PMTH2(:,:,:) = 4.*(MAX(ZZ_O_H,0.))**0.4*(ZZ_O_H-0.95)**2 END WHERE @@ -146,7 +149,16 @@ DO JK=IKTB+1,IKTE-1 END DO PMTH2(:,:,IKE)=PMTH2(:,:,IKE) * ZTSTAR(:,:)**2*ZWSTAR(:,:) PMTH2(:,:,KKU)=PMTH2(:,:,KKU) * ZTSTAR(:,:)**2*ZWSTAR(:,:) - +#else +WHERE(ZZ_O_H < 0.95 .AND. ZZ_O_H/=XUNDEF) + PMTH2(:,:,:) = 4.*(MAX(ZZ_O_H,0.))**0.4*BR_P2(ZZ_O_H-0.95) +END WHERE +DO JK=IKTB+1,IKTE-1 + PMTH2(:,:,JK) = PMTH2(:,:,JK) * BR_P2(ZTSTAR(:,:))*ZWSTAR(:,:) +END DO +PMTH2(:,:,IKE)=PMTH2(:,:,IKE) * BR_P2(ZTSTAR(:,:))*ZWSTAR(:,:) +PMTH2(:,:,KKU)=PMTH2(:,:,KKU) * BR_P2(ZTSTAR(:,:))*ZWSTAR(:,:) +#endif ! ! !* w'2th' @@ -155,11 +167,19 @@ PMWTH = 0. WHERE(ZZ_O_H <0.9 .AND. ZZ_O_H/=XUNDEF) PMWTH(:,:,:) = MAX(-7.9*(ABS(ZZ_O_H-0.35))**2.9 * (ABS(ZZ_O_H-1.))**0.58 + 0.37, 0.) END WHERE +#ifdef MNH_BITREP DO JK=IKTB+1,IKTE-1 PMWTH(:,:,JK) = PMWTH(:,:,JK) * ZWSTAR(:,:)**2*ZTSTAR(:,:) END DO PMWTH(:,:,IKE) = PMWTH(:,:,IKE) * ZWSTAR(:,:)**2*ZTSTAR(:,:) PMWTH(:,:,KKU) = PMWTH(:,:,KKU) * ZWSTAR(:,:)**2*ZTSTAR(:,:) +#else +DO JK=IKTB+1,IKTE-1 + PMWTH(:,:,JK) = PMWTH(:,:,JK) * BR_P2(ZWSTAR(:,:))*ZTSTAR(:,:) +END DO +PMWTH(:,:,IKE) = PMWTH(:,:,IKE) * BR_P2(ZWSTAR(:,:))*ZTSTAR(:,:) +PMWTH(:,:,KKU) = PMWTH(:,:,KKU) * BR_P2(ZWSTAR(:,:))*ZTSTAR(:,:) +#endif ! !---------------------------------------------------------------------------- END SUBROUTINE TM06 diff --git a/src/MNH/tridiag_thermo.f90 b/src/MNH/tridiag_thermo.f90 index 6f7cc20793cf0c2f89961783b1b210b5bb38f9bc..b1aedcdef303d59af453c48212951eb0ec485c19 100644 --- a/src/MNH/tridiag_thermo.f90 +++ b/src/MNH/tridiag_thermo.f90 @@ -160,6 +160,9 @@ USE MODI_SHUMAN #else USE MODI_SHUMAN_DEVICE #endif +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -219,7 +222,11 @@ ZMZM_RHODJ = MZM(KKA,KKU,KKL,PRHODJ) CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ) #endif !$acc kernels async +#ifdef MNH_BITREP ZRHODJ_DFDDTDZ_O_DZ2 = ZMZM_RHODJ*PDFDDTDZ/PDZZ**2 +#else +ZRHODJ_DFDDTDZ_O_DZ2 = ZMZM_RHODJ*PDFDDTDZ/BR_P2(PDZZ) +#endif !$acc end kernels ! !$acc kernels async diff --git a/src/MNH/tridiag_w.f90 b/src/MNH/tridiag_w.f90 index 8050ede969321fc4019f0ed2efe06ff65e715229..a518d51d51e5e2f7df4d82cd5c3d37c3b5bbffed 100644 --- a/src/MNH/tridiag_w.f90 +++ b/src/MNH/tridiag_w.f90 @@ -156,6 +156,9 @@ USE MODI_SHUMAN #else USE MODI_SHUMAN_DEVICE #endif +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -203,7 +206,11 @@ ZMZM_RHODJ = MZM(1,IKU,1,PRHODJ) CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ) #endif !$acc kernels async +#ifndef MNH_BITREP ZRHODJ_DFDDWDZ_O_DZ2 = PRHODJ*PDFDDWDZ/PMZF_DZZ**2 +#else +ZRHODJ_DFDDWDZ_O_DZ2 = PRHODJ*PDFDDWDZ/BR_P2(PMZF_DZZ) +#endif !$acc end kernels ! !$acc kernels async @@ -222,10 +229,17 @@ ZY=0. !! y(k) = (PRHODJ(k)+PRHODJ(k-1))/2.*PVARM(k)/PTSTEP !! - PRHODJ(k) * PF(k )/PMZF_PDZZ(k ) !! + PRHODJ(k-1) * PF(k-1)/PMZF_PDZZ(k-1) +!!#ifdef MNH_BITREP !! + PRHODJ(k) * PDFDDWDZ(k) * PVARM(k+1)/PMZF_DZZ(k)**2 !! - PRHODJ(k) * PDFDDWDZ(k) * PVARM(k) /PMZF_DZZ(k)**2 !! - PRHODJ(k-1) * PDFDDWDZ(k-1) * PVARM(k) /PMZF_DZZ(k-1)**2 !! + PRHODJ(k-1) * PDFDDWDZ(k-1) * PVARM(k-1)/PMZF_DZZ(k-1)**2 +!!#else +!! + PRHODJ(k) * PDFDDWDZ(k) * PVARM(k+1)/BR_P2(PMZF_DZZ(k)) +!! - PRHODJ(k) * PDFDDWDZ(k) * PVARM(k) /BR_P2(PMZF_DZZ(k)) +!! - PRHODJ(k-1) * PDFDDWDZ(k-1) * PVARM(k) /BR_P2(PMZF_DZZ(k-1)) +!! + PRHODJ(k-1) * PDFDDWDZ(k-1) * PVARM(k-1)/BR_P2(PMZF_DZZ(k-1)) +!!#endif ! !$acc kernels async ZY(:,:,IKB) = ZMZM_RHODJ(:,:,IKB)*PVARM(:,:,IKB)/PTSTEP & diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90 index b7973ae35387fe506912cbf4d30415e330d8c403..8c6185b6d23dfd6857b5801703d10a094d64983b 100644 --- a/src/MNH/turb.f90 +++ b/src/MNH/turb.f90 @@ -383,6 +383,7 @@ USE MODI_SECOND_MNH #ifdef MNH_BITREP USE MODI_BITREP #endif +USE MODE_MPPDB ! !! use, intrinsic :: ISO_C_BINDING ! @@ -877,7 +878,11 @@ IF (ORMC01) THEN PRINT *,'OPENACC: TURB::ORMC01 not yet implemented' CALL ABORT #endif +#ifndef MNH_BITREP ZUSTAR=(PSFU**2+PSFV**2)**(0.25) +#else + ZUSTAR=BR_POW(BR_P2(PSFU)+BR_P2(PSFV),0.25) +#endif IF (KRR>0) THEN ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV) ELSE @@ -924,9 +929,13 @@ END IF !* 4.2 compute the proportionality coefficient between wind and stress ! !$acc kernels - ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / & - (XMNH_TINY + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) & - ) +#ifdef MNH_BITREP + ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / & + (XMNH_TINY + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) ) +#else + ZCDUEFF(:,:) =-SQRT ( (BR_P2(PSFU(:,:)) + BR_P2(PSFV(:,:))) / & + (XMNH_TINY + BR_P2(ZUSLOPE(:,:)) + BR_P2(ZVSLOPE(:,:)) ) ) +#endif !$acc end kernels ! !* 4.6 compute the surface tangential fluxes @@ -1027,6 +1036,7 @@ CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI, & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS, & PDYP,PTHP,PSIGS,PWTH,PWRC,PWSV ) !$acc update self(PWTH,PWRC,PWSV) +CALL MPPDB_CHECK3DM("turb after turb_ver:PDYP,PTHP",PRECISION,PDYP,PTHP) ! #ifdef _OPENACC IF ( ( LBUDGET_TH .AND. ( ( KRRI >= 1 .AND. KRRL >= 1 ) .OR. ( KRRL >= 1 ) ) ) .OR. & @@ -1096,6 +1106,7 @@ IF (HTURBDIM=='3DIM') THEN PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS ) END IF !$acc update self(PSIGS,PRUS,PRVS,PRWS,PRSVS) +CALL MPPDB_CHECK3DM("turb after turb_hor_splt:PDYP,PTHP",PRECISION,PDYP,PTHP) ! ! #ifdef _OPENACC @@ -1163,6 +1174,7 @@ CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,PLEM,ZLEPS,PDYP,ZTRH, & ! !$acc update self(PTR,PDISS) !$acc update self(PDYP,PTHP) +CALL MPPDB_CHECK3DM("turb after tke_eps_sources:PDYP,PTHP",PRECISION,PDYP,PTHP) ! !$acc update self(PRTKES) IF (LBUDGET_TH) THEN @@ -1568,7 +1580,11 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT ( 1. + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) ) * & ( & ZRVSAT(:,:,:) * (1. + ZRVSAT(:,:,:)/ZEPS) & +#ifdef MNH_BITREP * ( -2.*PBETA/PT(:,:,:) + PGAM ) / PT(:,:,:)**2 & +#else + * ( -2.*PBETA/PT(:,:,:) + PGAM ) / BR_P2(PT(:,:,:)) & +#endif +ZDRVSATDT(:,:,:) * (1. + 2. * ZRVSAT(:,:,:)/ZEPS) & * ( PBETA/PT(:,:,:) - PGAM ) / PT(:,:,:) & ) & diff --git a/src/MNH/turb_hor_dyn_corr.f90 b/src/MNH/turb_hor_dyn_corr.f90 index 1a3a6188e1b170bc653998b4d6ca588abf44778b..0ec691e5f96325db6d4170882d6746fb7fb591c2 100644 --- a/src/MNH/turb_hor_dyn_corr.f90 +++ b/src/MNH/turb_hor_dyn_corr.f90 @@ -178,6 +178,9 @@ USE MODI_TRIDIAG_W ! USE MODI_SECOND_MNH USE MODE_MPPDB +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -307,7 +310,11 @@ IKU = SIZE(PUM,3) ! ! !$acc kernels async(1) +#ifndef MNH_BITREP ZDIRSINZW(:,:) = SQRT( 1. - PDIRCOSZW(:,:)**2 ) +#else +ZDIRSINZW(:,:) = SQRT( 1. - BR_P2(PDIRCOSZW(:,:)) ) +#endif !$acc end kernels ! #ifndef _OPENACC @@ -504,15 +511,25 @@ ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) & !$ acc wait(1) ! !$acc kernels async(4) -ZFLX(:,:,IKB-1) = & +#ifndef MNH_BITREP +ZFLX(:,:,IKB-1) = & PTAU11M(:,:) * PCOSSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2 & -2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) & + PTAU22M(:,:) * PSINSLOPE(:,:)**2 & + PTAU33M(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2 & +2. * PCDUEFF(:,:) * ( & PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) & - - PUSLOPEM(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) & - ) + - PUSLOPEM(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) ) +#else +ZFLX(:,:,IKB-1) = & + PTAU11M(:,:) * BR_P2(PCOSSLOPE(:,:)) * BR_P2(PDIRCOSZW(:,:)) & + -2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) & + + PTAU22M(:,:) * BR_P2(PSINSLOPE(:,:)) & + + PTAU33M(:,:) * BR_P2(PCOSSLOPE(:,:)) * BR_P2(ZDIRSINZW(:,:)) & + +2. * PCDUEFF(:,:) * ( & + PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) & + - PUSLOPEM(:,:) * BR_P2(PCOSSLOPE(:,:)) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) ) +#endif !$acc end kernels ! !!! wait for the computation of ZFLX(:,:,IKB) and ZFLX(:,:,IKB-1) @@ -695,15 +712,25 @@ ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) & ! ! extrapolates this flux under the ground with the surface flux !$acc kernels async(3) -ZFLX(:,:,IKB-1) = & +#ifndef MNH_BITREP +ZFLX(:,:,IKB-1) = & PTAU11M(:,:) * PSINSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2 & +2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) & + PTAU22M(:,:) * PCOSSLOPE(:,:)**2 & + PTAU33M(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2 & -2. * PCDUEFF(:,:)* ( & PUSLOPEM(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) & - + PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) & - ) + + PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) ) +#else +ZFLX(:,:,IKB-1) = & + PTAU11M(:,:) * BR_P2(PSINSLOPE(:,:)) * BR_P2(PDIRCOSZW(:,:)) & + +2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) & + + PTAU22M(:,:) * BR_P2(PCOSSLOPE(:,:)) & + + PTAU33M(:,:) * BR_P2(PSINSLOPE(:,:)) * BR_P2(ZDIRSINZW(:,:)) & + -2. * PCDUEFF(:,:)* ( & + PUSLOPEM(:,:) * BR_P2(PSINSLOPE(:,:)) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) & + + PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) ) +#endif !$acc end kernels ! !$acc kernels async(3) @@ -878,8 +905,13 @@ ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) & ! extrapolates this flux under the ground with the surface flux !$acc kernels async(3) ZFLX(:,:,IKB-1) = & - PTAU11M(:,:) * ZDIRSINZW(:,:)**2 & +#ifndef MNH_BITREP + PTAU11M(:,:) * ZDIRSINZW(:,:)**2 & + PTAU33M(:,:) * PDIRCOSZW(:,:)**2 & +#else + PTAU11M(:,:) * BR_P2(ZDIRSINZW(:,:)) & + + PTAU33M(:,:) * BR_P2(PDIRCOSZW(:,:)) & +#endif +2. * PCDUEFF(:,:)* PUSLOPEM(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) !$acc end kernels ! diff --git a/src/MNH/turb_hor_sv_corr.f90 b/src/MNH/turb_hor_sv_corr.f90 index 2f2f8f7b5cc0af628e82149a9ef4b7c68e45f9af..cc9d894a35bcf3b2138e524937f478b49b9202b2 100644 --- a/src/MNH/turb_hor_sv_corr.f90 +++ b/src/MNH/turb_hor_sv_corr.f90 @@ -106,6 +106,9 @@ USE MODI_EMOIST USE MODI_ETHETA ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -166,11 +169,20 @@ DO JSV=1,NSV IF (LLES_CALL) THEN IF (.NOT. L2D) THEN ZFLX(:,:,:) = XCHF / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP ( GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2 & + GY_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)**2 ) +#else + ( BR_P2(GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)) & + + BR_P2(GY_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)) ) +#endif ELSE ZFLX(:,:,:) = XCHF / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2 +#else + BR_P2(GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)) +#endif END IF CALL LES_MEAN_SUBGRID( -2.*ZCSVD*SQRT(PTKEM)*ZFLX/PLEPS, & X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV), .TRUE. ) diff --git a/src/MNH/turb_hor_thermo_corr.f90 b/src/MNH/turb_hor_thermo_corr.f90 index 2ac18385046c5ae6b39f1de50fbe111a475e2f66..66310b6477cd1223dbd456e94eacd33b737ba363 100644 --- a/src/MNH/turb_hor_thermo_corr.f90 +++ b/src/MNH/turb_hor_thermo_corr.f90 @@ -159,6 +159,9 @@ USE MODI_EMOIST USE MODI_ETHETA ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -281,22 +284,39 @@ IF ( ( KRRL > 0 .AND. OSUBG_COND) .OR. ( OTURB_FLX .AND. OCLOSE_OUT ) & #ifndef _OPENACC IF (.NOT. L2D) THEN ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP ( GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)**2 + GY_M_M(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY)**2 ) +#else + ( BR_P2(GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)) + & + BR_P2(GY_M_M(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY)) ) +#endif ELSE ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)**2 +#else + BR_P2(GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)) +#endif END IF #else IF (.NOT. L2D) THEN CALL GX_M_M_DEVICE(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) CALL GY_M_M_DEVICE(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY,ZTMP2_DEVICE) !$acc kernels +#ifdef MNH_BITREP ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * ( ZTMP1_DEVICE**2 + ZTMP2_DEVICE**2 ) +#else + ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * ( BR_P2(ZTMP1_DEVICE) + BR_P2(ZTMP2_DEVICE) ) +#endif !$acc end kernels ELSE CALL GX_M_M_DEVICE(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE) !$acc kernels +#ifdef MNH_BITREP ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * ZTMP1_DEVICE**2 +#else + ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * BR_P2(ZTMP1_DEVICE) +#endif !$acc end kernels END IF #endif @@ -627,23 +647,40 @@ IF ( ( KRRL > 0 .AND. OSUBG_COND) .OR. ( OTURB_FLX .AND. OCLOSE_OUT ) & #ifndef _OPENACC IF (.NOT. L2D) THEN ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP ( GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 + & GY_M_M(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY)**2 ) +#else + ( BR_P2(GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)) + & + BR_P2(GY_M_M(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY)) ) +#endif ELSE ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * & +#ifdef MNH_BITREP ( GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 ) +#else + ( BR_P2(GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)) ) +#endif END IF #else IF (.NOT. L2D) THEN CALL GX_M_M_DEVICE(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE) CALL GY_M_M_DEVICE(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZTMP2_DEVICE) !$acc kernels +#ifdef MNH_BITREP ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * ( ZTMP1_DEVICE**2 + ZTMP2_DEVICE**2 ) +#else + ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * ( BR_P2(ZTMP1_DEVICE) + BR_P2(ZTMP2_DEVICE) ) +#endif !$acc end kernels ELSE CALL GX_M_M_DEVICE(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE) !$acc kernels +#ifdef MNH_BITREP ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * ZTMP1_DEVICE**2 +#else + ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * BR_P2(ZTMP1_DEVICE) +#endif !$acc end kernels END IF #endif diff --git a/src/MNH/turb_hor_uv.f90 b/src/MNH/turb_hor_uv.f90 index bb83050657fd1a4e6e391044ea4063ba1a0139ba..28bdeaeb8bff8cc8766d9fe587cf53b8c851bc2d 100644 --- a/src/MNH/turb_hor_uv.f90 +++ b/src/MNH/turb_hor_uv.f90 @@ -158,6 +158,9 @@ USE MODI_COEFJ USE MODI_LES_MEAN_SUBGRID ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -264,7 +267,11 @@ IKU = SIZE(PUM,3) ILENCH=LEN(YCOMMENT) ! !$acc kernels +#ifdef MNH_BITREP ZDIRSINZW(:,:) = SQRT( 1. - PDIRCOSZW(:,:)**2 ) +#else +ZDIRSINZW(:,:) = SQRT( 1. - BR_P2(PDIRCOSZW(:,:)) ) +#endif !$acc end kernels ! #ifndef _OPENACC @@ -376,6 +383,7 @@ ZFLX(:,:,IKB) = - XCMFS * ZTMP2_DEVICE(:,:,1) * ( ZTMP5_DEVICE(:,:,1) + ZTMP6_DE ! ! extrapolates this flux under the ground with the surface flux !$acc kernels +#ifdef MNH_BITREP ZFLX(:,:,IKB-1) = & PTAU11M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:)**2 & +PTAU12M(:,:) * (PCOSSLOPE(:,:)**2 - PSINSLOPE(:,:)**2) * & @@ -385,8 +393,19 @@ ZFLX(:,:,IKB-1) = & -PCDUEFF(:,:) * ( & 2. * PUSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * & PDIRCOSZW(:,:) * ZDIRSINZW(:,:) & - +PVSLOPEM(:,:) * (PCOSSLOPE(:,:)**2 - PSINSLOPE(:,:)**2) * ZDIRSINZW(:,:) & - ) + +PVSLOPEM(:,:) * (PCOSSLOPE(:,:)**2 - PSINSLOPE(:,:)**2) * ZDIRSINZW(:,:) ) +#else +ZFLX(:,:,IKB-1) = & + PTAU11M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * BR_P2(PDIRCOSZW(:,:)) & + +PTAU12M(:,:) * (BR_P2(PCOSSLOPE(:,:)) - BR_P2(PSINSLOPE(:,:))) * & + BR_P2(PDIRCOSZW(:,:)) & + -PTAU22M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) & + +PTAU33M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * BR_P2(ZDIRSINZW(:,:)) & + -PCDUEFF(:,:) * ( & + 2. * PUSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * & + PDIRCOSZW(:,:) * ZDIRSINZW(:,:) & + +PVSLOPEM(:,:) * (BR_P2(PCOSSLOPE(:,:)) - BR_P2(PSINSLOPE(:,:))) * ZDIRSINZW(:,:) ) +#endif !$acc end kernels ! #ifndef _OPENACC diff --git a/src/MNH/turb_ver_dyn_flux.f90 b/src/MNH/turb_ver_dyn_flux.f90 index 0304e7e5646154dfb85c3d7fbf2566b55dc71ed5..04df9fa5daf05967991a67e875f2bbf060943a23 100644 --- a/src/MNH/turb_ver_dyn_flux.f90 +++ b/src/MNH/turb_ver_dyn_flux.f90 @@ -322,6 +322,9 @@ USE MODI_LES_MEAN_SUBGRID ! USE MODI_SECOND_MNH USE MODE_ll +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -454,7 +457,11 @@ IKTE=IKT-JPVEXT_TURB !$acc kernels ZSOURCE(:,:,:) = 0. ! +#ifndef MNH_BITREP ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2) +#else +ZDIRSINZW(:,:) = SQRT(1.-BR_P2(PDIRCOSZW(:,:))) +#endif ! compute the coefficients for the uncentred gradient computation near the ! ground ! @@ -484,14 +491,22 @@ ZVSLOPEM(:,:,1)=PVSLOPEM(:,:) #ifndef _OPENACC ZA(:,:,:) = -PTSTEP * XCMFS * & MXM( ZKEFF ) * MXM(MZM(KKA,KKU,KKL, PRHODJ )) / & +#ifndef MNH_BITREP MXM( PDZZ )**2 +#else + BR_P2(MXM( PDZZ )) +#endif #else CALL MXM_DEVICE( ZKEFF, ZTMP1_DEVICE ) CALL MZM_DEVICE(PRHODJ, ZTMP2_DEVICE ) CALL MXM_DEVICE(ZTMP2_DEVICE,ZTMP3_DEVICE) CALL MXM_DEVICE( PDZZ, ZTMP4_DEVICE ) !$acc kernels +#ifndef MNH_BITREP ZA(:,:,:) = -PTSTEP * XCMFS * ZTMP1_DEVICE * ZTMP3_DEVICE / ZTMP4_DEVICE**2 +#else +ZA(:,:,:) = -PTSTEP * XCMFS * ZTMP1_DEVICE * ZTMP3_DEVICE / BR_P2(ZTMP4_DEVICE) +#endif !$acc end kernels #endif ! @@ -506,7 +521,11 @@ ENDIF ! compute the coefficient between the vertical flux and the 2 components of the ! wind following the slope !$acc kernels +#ifndef MNH_BITREP ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) & +#else +ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (BR_P2(PDIRCOSZW(:,:)) - BR_P2(ZDIRSINZW(:,:))) & +#endif * PCOSSLOPE(:,:) ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PSINSLOPE(:,:) @@ -890,14 +909,22 @@ END IF #ifndef _OPENACC ZA(:,:,:) = - PTSTEP * XCMFS * & MYM( ZKEFF ) * MYM(MZM(KKA,KKU,KKL, PRHODJ )) / & +#ifndef MNH_BITREP MYM( PDZZ )**2 +#else + BR_P2(MYM( PDZZ )) +#endif #else CALL MYM_DEVICE( ZKEFF, ZTMP1_DEVICE ) CALL MYM_DEVICE( PDZZ, ZTMP2_DEVICE ) CALL MZM_DEVICE(PRHODJ, ZTMP3_DEVICE ) CALL MYM_DEVICE(ZTMP3_DEVICE,ZTMP4_DEVICE) !$acc kernels +#ifndef MNH_BITREP ZA(:,:,:) = - PTSTEP * XCMFS * ZTMP1_DEVICE * ZTMP4_DEVICE / ZTMP2_DEVICE**2 +#else +ZA(:,:,:) = - PTSTEP * XCMFS * ZTMP1_DEVICE * ZTMP4_DEVICE / BR_P2(ZTMP2_DEVICE) +#endif #endif ! ! @@ -906,7 +933,11 @@ IF(CPROGRAM/='AROME ') ZA(:,1,:)=ZA(:,IJE,:) ! Compute the source of V wind component ! compute the coefficient between the vertical flux and the 2 components of the ! wind following the slope +#ifndef MNH_BITREP ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) & +#else +ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (BR_P2(PDIRCOSZW(:,:)) - BR_P2(ZDIRSINZW(:,:))) & +#endif * PSINSLOPE(:,:) ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PCOSSLOPE(:,:) diff --git a/src/MNH/turb_ver_sv_corr.f90 b/src/MNH/turb_ver_sv_corr.f90 index 48c502605082e2c9fc700ff5c012d232c9b6e776..a8b026339f9aaae01a4cb57448f5ce7fa02272d8 100644 --- a/src/MNH/turb_ver_sv_corr.f90 +++ b/src/MNH/turb_ver_sv_corr.f90 @@ -120,6 +120,9 @@ USE MODI_ETHETA USE MODI_LES_MEAN_SUBGRID ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -183,7 +186,11 @@ DO JSV=1,NSV ! IF (LLES_CALL) THEN ! approximation: diagnosed explicitely (without implicit term) +#ifdef MNH_BITREP ZFLXZ(:,:,:) = PPSI_SV(:,:,:,JSV)*GZ_M_W(KKA,KKU,KKL,PSVM(:,:,:,JSV),PDZZ)**2 +#else + ZFLXZ(:,:,:) = PPSI_SV(:,:,:,JSV)*BR_P2(GZ_M_W(KKA,KKU,KKL,PSVM(:,:,:,JSV),PDZZ)) +#endif ZFLXZ(:,:,:) = XCHF / ZCSVD * PLM * PLEPS * MZF(KKA,KKU,KKL,ZFLXZ(:,:,:) ) CALL LES_MEAN_SUBGRID( -2.*ZCSVD*SQRT(PTKEM)*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV) ) CALL LES_MEAN_SUBGRID( MZF(KKA,KKU,KKL,PWM)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) ) diff --git a/src/MNH/turb_ver_sv_flux.f90 b/src/MNH/turb_ver_sv_flux.f90 index e11beff64a0cac292d0f24290637c2a392f489f1..068ba05f28483b6a5e5c79829f46993f01557521 100644 --- a/src/MNH/turb_ver_sv_flux.f90 +++ b/src/MNH/turb_ver_sv_flux.f90 @@ -286,6 +286,9 @@ USE MODE_FMWRIT USE MODI_LES_MEAN_SUBGRID ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -393,7 +396,11 @@ DO JSV=1,ISV ! Preparation of the arguments for TRIDIAG ZA(:,:,:) = -PTSTEP*XCHF*PPSI_SV(:,:,:,JSV) * & ZKEFF * MZM(KKA,KKU,KKL,PRHODJ) / & +#ifdef MNH_BITREP PDZZ**2 +#else + BR_P2(PDZZ) +#endif ! ! Compute the sources for the JSVth scalar variable diff --git a/src/MNH/turb_ver_thermo_corr.f90 b/src/MNH/turb_ver_thermo_corr.f90 index a54aae89a5a7907bd29c515c3d6aabcc6d449874..89b2dbe05cd7d099a0525ffdc8ab6eb2fa0a0aa3 100644 --- a/src/MNH/turb_ver_thermo_corr.f90 +++ b/src/MNH/turb_ver_thermo_corr.f90 @@ -350,6 +350,9 @@ USE MODI_LES_MEAN_SUBGRID USE MODE_PRANDTL ! USE MODI_SECOND_MNH +#ifdef MNH_BITREP +USE MODI_BITREP +#endif ! IMPLICIT NONE ! @@ -545,11 +548,19 @@ END IF ! ! Compute the turbulent variance F and F' at time t-dt. #ifndef _OPENACC +#ifdef MNH_BITREP ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(KKA,KKU,KKL,PPHI3*PDTH_DZ**2) +#else + ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(KKA,KKU,KKL,PPHI3*BR_P2(PDTH_DZ)) +#endif ZDFDDTDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately #else !$acc kernels +#ifdef MNH_BITREP ZTMP1_DEVICE = PPHI3*PDTH_DZ**2 +#else + ZTMP1_DEVICE = PPHI3*BR_P2(PDTH_DZ) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels @@ -697,6 +708,7 @@ END IF #endif ! ! special case near the ground ( uncentred gradient ) +#ifdef MNH_BITREP ZFLXZ(:,:,IKB) = XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB) & * PLEPS(:,:,IKB) & *( PEXPL * & @@ -708,13 +720,30 @@ END IF +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) & +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 & ) +#else + ZFLXZ(:,:,IKB) = XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB) & + * PLEPS(:,:,IKB) & + *( PEXPL * & + BR_P2( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) & + +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) & + +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) ) & + +PIMPL * & + BR_P2( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) & + +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) & + +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) ) & + ) +#endif ! ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB) ! ZFLXZ = MAX(0., ZFLXZ) ! IF (KRRL > 0) THEN +#ifdef MNH_BITREP PSIGS(:,:,:) = ZFLXZ(:,:,:) * PATHETA(:,:,:)**2 +#else + PSIGS(:,:,:) = ZFLXZ(:,:,:) * BR_P2(PATHETA(:,:,:)) +#endif END IF !$acc end kernels ! @@ -1101,10 +1130,18 @@ END IF ! ! Compute the turbulent variance F and F' at time t-dt. #ifndef _OPENACC +#ifdef MNH_BITREP ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(KKA,KKU,KKL,PPSI3*PDR_DZ**2) +#else + ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(KKA,KKU,KKL,PPSI3*BR_P2(PDR_DZ)) +#endif #else !$acc kernels +#ifdef MNH_BITREP ZTMP1_DEVICE = PPSI3*PDR_DZ**2 +#else + ZTMP1_DEVICE = PPSI3*BR_P2(PDR_DZ) +#endif !$acc end kernels CALL MZF_DEVICE(KKA,KKU,KKL,ZTMP1_DEVICE,ZTMP2_DEVICE) !$acc kernels @@ -1262,6 +1299,7 @@ END IF #endif ! ! special case near the ground ( uncentred gradient ) +#ifdef MNH_BITREP ZFLXZ(:,:,IKB) = XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB) & * PLEPS(:,:,IKB) & *( PEXPL * & @@ -1273,11 +1311,28 @@ END IF +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) & +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 & ) +#else + ZFLXZ(:,:,IKB) = XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB) & + * PLEPS(:,:,IKB) & + *( PEXPL * & + BR_P2( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) & + +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) & + +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) & + +PIMPL * & + BR_P2( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) & + +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) & + +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) & + ) +#endif ! ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB) ! IF ( KRRL > 0 ) THEN +#ifdef MNH_BITREP PSIGS(:,:,:) = PSIGS(:,:,:) + PAMOIST(:,:,:) **2 * ZFLXZ(:,:,:) +#else + PSIGS(:,:,:) = PSIGS(:,:,:) + BR_P2(PAMOIST(:,:,:)) * ZFLXZ(:,:,:) +#endif END IF !$acc end kernels ! stores <Rnp Rnp>