Ignore:
Timestamp:
Jun 17, 2025, 4:57:26 PM (7 days ago)
Author:
yann meurdesoif
Message:

GPU porting : revert of reverted commit from rev5511 (initial commit) and rev5561 (reverted commit).
This commit imply a lost of convergence in production mode due to arithmetics rounding effects (possibly by changing order of operations)

YM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90

    r5561 r5715  
    55
    66  USE clesphys_mod_h
    7     implicit none
     7  IMPLICIT NONE
     8  LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
     9  LOGICAL, SAVE :: firstcall = .TRUE.
     10  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
     11 
     12  INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
     13  INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
     14
    815
    916contains
    10 
    11   SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
    12        zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress)
    13 
    14     ! Parametrization of the momentum flux deposition due to a discrete
    15     ! number of gravity waves.
    16     ! Author: F. Lott, A. de la Camara
    17     ! July, 24th, 2014
    18     ! Gaussian distribution of the source, source is vorticity squared
    19     ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
    20     ! Lott et al (JAS, 2010, vol 67, page 157-170)
    21     ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
    22 
    23 !  ONLINE:
    24 USE yoegwd_mod_h
    25         USE yomcst_mod_h
    26 use dimphy, only: klon, klev
    27     use assert_m, only: assert
    28     USE ioipsl_getin_p_mod, ONLY : getin_p
    29     USE vertical_layers_mod, ONLY : presnivs
    30 
    31 
    32 !  OFFLINE:
    33 !   include "dimensions_mod.f90"
    34 !   include "dimphy.h"
    35 !END DIFFERENCE
    36 
    37     ! 0. DECLARATIONS:
    38 
    39     ! 0.1 INPUTS
    40     REAL, intent(in)::DTIME ! Time step of the Physics
    41     REAL, intent(in):: PP(:, :) ! (KLON, KLEV) Pressure at full levels
    42     REAL, intent(in):: ROT(:,:) ! Relative vorticity             
    43     REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
    44     REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
    45     REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
    46     REAL, intent(in):: PLAT(:) ! (KLON) LATITUDE                   
    47 
    48     ! 0.2 OUTPUTS
    49     REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
    50 
    51     REAL, intent(inout):: d_u(:, :), d_v(:, :)
    52     REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
    53     REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress
    54     ! (KLON, KLEV) tendencies on winds
    55 
    56     ! O.3 INTERNAL ARRAYS
    57     REAL BVLOW(klon)  !  LOW LEVEL BV FREQUENCY
    58     REAL ROTBA(KLON),CORIO(KLON)  !  BAROTROPIC REL. VORTICITY AND PLANETARY
    59     REAL UZ(KLON, KLEV + 1)
    60 
    61     INTEGER II, JJ, LL
    62 
    63     ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
    64 
    65     REAL DELTAT
    66 
    67     ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
    68 
    69     INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
    70     INTEGER JK, JP, JO, JW
    71     INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
    72     REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
    73     REAL CMIN, CMAX ! Min and Max absolute ph. vel.
    74     REAL CPHA ! absolute PHASE VELOCITY frequency
    75     REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
    76     REAL ZP(NW, KLON) ! Horizontal wavenumber angle
    77     REAL ZO(NW, KLON) ! Absolute frequency !
    78 
    79     ! Waves Intr. freq. at the 1/2 lev surrounding the full level
    80     REAL ZOM(NW, KLON), ZOP(NW, KLON)
    81 
    82     ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    83     REAL WWM(NW, KLON), WWP(NW, KLON)
    84 
    85     REAL RUW0(NW, KLON) ! Fluxes at launching level
    86 
    87     REAL RUWP(NW, KLON), RVWP(NW, KLON)
    88     ! Fluxes X and Y for each waves at 1/2 Levels
    89 
    90     INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
    91 
    92     REAL XLAUNCH ! Controle the launching altitude
    93     REAL XTROP ! SORT of Tropopause altitude
    94     REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
    95     REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
    96 
    97     REAL PRMAX ! Maximum value of PREC, and for which our linear formula
    98     ! for GWs parameterisation apply
    99 
    100     ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
    101 
    102     REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
    103     REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS
    104     REAL RUWFRT,SATFRT
    105 
    106     ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
    107 
    108     REAL H0 ! Characteristic Height of the atmosphere
    109     REAL DZ ! Characteristic depth of the source!
    110     REAL PR, TR ! Reference Pressure and Temperature
    111 
    112     REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
    113 
    114     REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
    115     REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
    116     REAL PSEC ! Security to avoid division by 0 pressure
    117     REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
    118     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
    119     REAL BVSEC ! Security to avoid negative BVF
    120 
    121     REAL, DIMENSION(klev+1) ::HREF
    122     LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
    123     LOGICAL, SAVE :: firstcall = .TRUE.
    124   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
    125 
    126     CHARACTER (LEN=20) :: modname='acama_gwd_rando_m'
     17 
     18  SUBROUTINE ACAMA_GWD_rando_first
     19  use dimphy, only: klev
     20  USE ioipsl_getin_p_mod, ONLY : getin_p
     21  IMPLICIT NONE
     22    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
    12723    CHARACTER (LEN=80) :: abort_message
    128 
    129 
    130 
    131   IF (firstcall) THEN
     24 
     25    IF (firstcall) THEN
    13226    ! Cle introduite pour resoudre un probleme de non reproductibilite
    13327    ! Le but est de pouvoir tester de revenir a la version precedenete
    … …  
    14034    firstcall=.false.
    14135!    CALL iophys_ini(dtime)
    142   ENDIF
     36    ENDIF
     37  END SUBROUTINE ACAMA_GWD_rando_first
     38
     39  SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
     40       zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress)
     41
     42    ! Parametrization of the momentum flux deposition due to a discrete
     43    ! number of gravity waves.
     44    ! Author: F. Lott, A. de la Camara
     45    ! July, 24th, 2014
     46    ! Gaussian distribution of the source, source is vorticity squared
     47    ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
     48    ! Lott et al (JAS, 2010, vol 67, page 157-170)
     49    ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
     50
     51!  ONLINE:
     52USE yoegwd_mod_h
     53        USE yomcst_mod_h
     54use dimphy, only: klon, klev
     55    use assert_m, only: assert
     56    USE vertical_layers_mod, ONLY : presnivs
     57    USE lmdz_fake_call, ONLY : fake_call
     58
     59!  OFFLINE:
     60!   include "dimensions_mod.f90"
     61!   include "dimphy.h"
     62!END DIFFERENCE
     63
     64    ! 0. DECLARATIONS:
     65
     66    ! 0.1 INPUTS
     67    REAL, intent(in)::DTIME ! Time step of the Physics
     68    REAL, intent(in):: PP(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels
     69    REAL, intent(in):: ROT(KLON,KLEV) ! Relative vorticity             
     70    REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels
     71    REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels
     72    REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels
     73    REAL, intent(in):: PLAT(KLON) ! (KLON) LATITUDE                   
     74
     75    ! 0.2 OUTPUTS
     76    REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses
     77
     78    REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV)
     79    REAL, intent(inout):: east_gwstress(KLON, KLEV) !  Profile of eastward stress
     80    REAL, intent(inout):: west_gwstress(KLON, KLEV) !  Profile of westward stress
     81    ! (KLON, KLEV) tendencies on winds
     82
     83    ! O.3 INTERNAL ARRAYS
     84    REAL BVLOW(klon)  !  LOW LEVEL BV FREQUENCY
     85    REAL ROTBA(KLON),CORIO(KLON)  !  BAROTROPIC REL. VORTICITY AND PLANETARY
     86    REAL UZ(KLON, KLEV + 1)
     87
     88    INTEGER II, JJ, LL
     89
     90    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
     91
     92    REAL DELTAT
     93
     94    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
     95
     96    INTEGER JK, JP, JO, JW
     97    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
     98    REAL CMIN, CMAX ! Min and Max absolute ph. vel.
     99    REAL CPHA ! absolute PHASE VELOCITY frequency
     100    REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude
     101    REAL ZP(KLON, NW) ! Horizontal wavenumber angle
     102    REAL ZO(KLON,NW) ! Absolute frequency !
     103
     104    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
     105    REAL ZOM(KLON, NW), ZOP(KLON, NW)
     106
     107    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
     108    REAL WWM(KLON, NW), WWP(KLON, NW)
     109
     110    REAL RUW0(KLON, NW) ! Fluxes at launching level
     111
     112    REAL RUWP(KLON, NW), RVWP(KLON, NW)
     113    ! Fluxes X and Y for each waves at 1/2 Levels
     114
     115    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
     116
     117    REAL XLAUNCH ! Controle the launching altitude
     118    REAL XTROP ! SORT of Tropopause altitude
     119    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
     120    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
     121
     122    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
     123    ! for GWs parameterisation apply
     124
     125    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
     126
     127    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
     128    REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS
     129    REAL RUWFRT,SATFRT
     130
     131    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
     132
     133    REAL H0 ! Characteristic Height of the atmosphere
     134    REAL DZ ! Characteristic depth of the source!
     135    REAL PR, TR ! Reference Pressure and Temperature
     136
     137    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
     138
     139    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
     140    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
     141    REAL PSEC ! Security to avoid division by 0 pressure
     142    REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
     143    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
     144    REAL BVSEC ! Security to avoid negative BVF
     145
     146    REAL, DIMENSION(klev+1) ::HREF
     147
     148    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
     149    CHARACTER (LEN=80) :: abort_message
     150
    143151
    144152    !-----------------------------------------------------------------
    … …  
    211219!  END ONLINE
    212220
     221    CALL FAKE_CALL(BVLOW) ! to be suppress in future
     222    CALL FAKE_CALL(CORIO) ! to be suppress in future
     223    CALL FAKE_CALL(ROTBA) ! to be suppress in future
     224   
    213225    IF(DELTAT < DTIME)THEN
    214226!       PRINT *, 'flott_gwd_rando: deltat < dtime!'
    … …  
    276288            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
    277289    end DO
    278     BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
     290    BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
    279291         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
    280292         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
    … …  
    345357             DO II = 1, KLON
    346358                ! Angle (0 or PI so far)
    347                 ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
     359                ! ZP(II,JW) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
    348360                !      * RPI / 2.
    349361                ! Angle between 0 and pi
    350                   ZP(JW, II) = MOD(TT(II, JW) * 10., 1.) * RPI
     362                  ZP(II,JW) = MOD(TT(II, JW) * 10., 1.) * RPI
    351363! TEST WITH POSITIVE WAVES ONLY (Part I/II)
    352 !               ZP(JW, II) = 0.
     364!               ZP(II,JW) = 0.
    353365                ! Horizontal wavenumber amplitude
    354                 ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
     366                ZK(II,JW) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
    355367                ! Horizontal phase speed
    356368                CPHA = 0.
    … …  
    361373                IF (CPHA.LT.0.)  THEN
    362374                   CPHA = -1.*CPHA
    363                    ZP(JW,II) = ZP(JW,II) + RPI
     375                   ZP(II,JW) = ZP(II,JW) + RPI
    364376! TEST WITH POSITIVE WAVES ONLY (Part II/II)
    365 !               ZP(JW, II) = 0.
     377!               ZP(II,JW) = 0.
    366378                ENDIF
    367379                CPHA = CPHA + CMIN !we dont allow |c|<1m/s
    368380                ! Absolute frequency is imposed
    369                 ZO(JW, II) = CPHA * ZK(JW, II)
     381                ZO(II, JW) = CPHA * ZK(II,JW)
    370382                ! Intrinsic frequency is imposed
    371                 ZO(JW, II) = ZO(JW, II) &
    372                      + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
    373                      + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
     383                ZO(II, JW) = ZO(II, JW) &
     384                     + ZK(II,JW) * COS(ZP(II,JW)) * UH(II, LAUNCH) &
     385                     + ZK(II,JW) * SIN(ZP(II,JW)) * VH(II, LAUNCH)
    374386                ! Momentum flux at launch lev
    375387                ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE
    376388                !  RIGHT IN THE SH (GWD4 after 1990)
    377                   RUW0(JW, II) = 0.
     389                  RUW0(II, JW) = 0.
    378390                 DO JJ = 1, NA
    379                     RUW0(JW, II) = RUW0(JW,II) + &
     391                    RUW0(II, JW) = RUW0(II, JW) + &
    380392         2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.)
    381393                END DO
    382                 RUW0(JW, II) = RUWFRT &
    383                           * EXP(RUW0(JW,II))/1250. &  ! 2 mpa at south pole
     394                RUW0(II,JW) = RUWFRT &
     395                          * EXP(RUW0(II,JW))/1250. &  ! 2 mpa at south pole
    384396       *((1.05+SIN(PLAT(II)*RPI/180.))/(1.01+SIN(PLAT(II)*RPI/180.))-2.05/2.01)
    385                 ! RUW0(JW, II) = RUWFRT
     397                ! RUW0(II,JW) = RUWFRT
    386398             ENDDO
    387399    end DO
    … …  
    397409
    398410       ! Evaluate intrinsic frequency at launching altitude:
    399        ZOP(JW, :) = ZO(JW, :) &
    400             - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
    401             - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
     411       ZOP(:, JW) = ZO(:, JW) &
     412            - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LAUNCH) &
     413            - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LAUNCH)
    402414
    403415       ! VERSION WITH FRONTAL SOURCES
    … …  
    406418
    407419       ! tanh limitation for values above CORIO (inertial instability).
    408        ! WWP(JW, :) = RUW0(JW, :) &
    409        WWP(JW, :) = RUWFRT      &
     420       ! WWP(:,JW) = RUW0(:,JW) &
     421       WWP(:,JW) = RUWFRT      &
    410422       !     * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 &
    411423       !    * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) &
    … …  
    413425       !    * (CORIO(:)*CORIO(:)) &
    414426       ! MODERATION BY THE DEPTH OF THE SOURCE (DZ HERE)
    415        !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(JW, :)),ZOISEC)**2 &
    416        !      *ZK(JW, :)**2*DZ**2) &
     427       !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(:,JW)),ZOISEC)**2 &
     428       !      *ZK(:,JW)**2*DZ**2) &
    417429       ! COMPLETE FORMULA:
    418430            !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) &
    … …  
    420432       !  RESTORE DIMENSION OF A FLUX
    421433       !     *RD*TR/PR
    422        !     *1. + RUW0(JW, :)
     434       !     *1. + RUW0(:,JW)
    423435             *1.
    424436
    … …  
    429441       ! Put the stress in the right direction:
    430442
    431         RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    432         RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     443        RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     444        RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    433445
    434446    end DO
    … …  
    440452       RVW(:, LL) = 0
    441453       DO JW = 1, NW
    442           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    443           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
     454          RUW(:, LL) = RUW(:, LL) + RUWP(:,JW)
     455          RVW(:, LL) = RVW(:, LL) + RVWP(:,JW)
    444456       end DO
    445457    end DO
    … …  
    453465       ! to the next)
    454466       DO JW = 1, NW
    455           ZOM(JW, :) = ZOP(JW, :)
    456           WWM(JW, :) = WWP(JW, :)
     467          ZOM(:, JW) = ZOP(:,JW)
     468          WWM(:,JW) = WWP(:,JW)
    457469          ! Intrinsic Frequency
    458           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
    459                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
     470          ZOP(:,JW) = ZO(:, JW) - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LL + 1) &
     471               - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LL + 1)
    460472
    461473          ! No breaking (Eq.6)
    462474          ! Dissipation (Eq. 8)
    463           WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
     475          WWP(:,JW) = WWM(:,JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
    464476               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
    465                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
    466                * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
     477               / MAX(ABS(ZOP(:,JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
     478               * ZK(:,JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
    467479
    468480          ! Critical levels (forced to zero if intrinsic frequency changes sign)
    469481          ! Saturation (Eq. 12)
    470           WWP(JW, :) = min(WWP(JW, :), MAX(0., &
    471                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
     482          WWP(:,JW) = min(WWP(:,JW), MAX(0., &
     483               SIGN(1., ZOP(:,JW) * ZOM(:, JW))) * ABS(ZOP(:,JW))**3 &
    472484          !    / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 &
    473485               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 &
    474486!              *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 &
    475487               *SATFRT**2       &
    476                / ZK(JW, :)**4)
     488               / ZK(:,JW)**4)
    477489       end DO
    478490
    … …  
    481493
    482494       DO JW = 1, NW
    483           RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    484           RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     495          RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     496          RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    485497       end DO
    486498
    … …  
    489501
    490502       DO JW = 1, NW
    491           RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
    492           RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
    493           EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW)
    494           WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW)
     503          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:,JW)
     504          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:,JW)
     505          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:,JW))/REAL(NW)
     506          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:,JW))/REAL(NW)
    495507       end DO
    496508    end DO
Note: See TracChangeset for help on using the changeset viewer.