Ignore:
Timestamp:
Jul 28, 2025, 7:23:15 PM (6 days ago)
Author:
aborella
Message:

Merge with trunk r5789

Location:
LMDZ6/branches/contrails
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails

  • LMDZ6/branches/contrails/libf/phylmd/acama_gwd_rando_m.f90

    r5309 r5791  
    22! $Id$
    33!
     4!$gpum horizontal klon
    45module ACAMA_GWD_rando_m
    56
    67  USE clesphys_mod_h
    7     implicit none
     8  IMPLICIT NONE
     9  LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
     10  LOGICAL, SAVE :: firstcall = .TRUE.
     11  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
     12 
     13  INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
     14  INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
     15
    816
    917contains
    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'
     18 
     19  SUBROUTINE ACAMA_GWD_rando_first
     20  use dimphy, only: klev
     21  USE ioipsl_getin_p_mod, ONLY : getin_p
     22  IMPLICIT NONE
     23    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
    12724    CHARACTER (LEN=80) :: abort_message
    128 
    129 
    130 
    131   IF (firstcall) THEN
     25 
     26    IF (firstcall) THEN
    13227    ! Cle introduite pour resoudre un probleme de non reproductibilite
    13328    ! Le but est de pouvoir tester de revenir a la version precedenete
    … …  
    14035    firstcall=.false.
    14136!    CALL iophys_ini(dtime)
    142   ENDIF
     37    ENDIF
     38  END SUBROUTINE ACAMA_GWD_rando_first
     39
     40  SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
     41       zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress)
     42
     43    ! Parametrization of the momentum flux deposition due to a discrete
     44    ! number of gravity waves.
     45    ! Author: F. Lott, A. de la Camara
     46    ! July, 24th, 2014
     47    ! Gaussian distribution of the source, source is vorticity squared
     48    ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
     49    ! Lott et al (JAS, 2010, vol 67, page 157-170)
     50    ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
     51
     52!  ONLINE:
     53USE yoegwd_mod_h
     54        USE yomcst_mod_h
     55use dimphy, only: klon, klev
     56    use assert_m, only: assert
     57    USE vertical_layers_mod, ONLY : presnivs
     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    !-----------------------------------------------------------------
    … …  
    276284            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
    277285    end DO
    278     BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
     286    BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
    279287         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
    280288         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
    … …  
    345353             DO II = 1, KLON
    346354                ! Angle (0 or PI so far)
    347                 ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
     355                ! ZP(II,JW) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
    348356                !      * RPI / 2.
    349357                ! Angle between 0 and pi
    350                   ZP(JW, II) = MOD(TT(II, JW) * 10., 1.) * RPI
     358                  ZP(II,JW) = MOD(TT(II, JW) * 10., 1.) * RPI
    351359! TEST WITH POSITIVE WAVES ONLY (Part I/II)
    352 !               ZP(JW, II) = 0.
     360!               ZP(II,JW) = 0.
    353361                ! Horizontal wavenumber amplitude
    354                 ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
     362                ZK(II,JW) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
    355363                ! Horizontal phase speed
    356364                CPHA = 0.
    … …  
    361369                IF (CPHA.LT.0.)  THEN
    362370                   CPHA = -1.*CPHA
    363                    ZP(JW,II) = ZP(JW,II) + RPI
     371                   ZP(II,JW) = ZP(II,JW) + RPI
    364372! TEST WITH POSITIVE WAVES ONLY (Part II/II)
    365 !               ZP(JW, II) = 0.
     373!               ZP(II,JW) = 0.
    366374                ENDIF
    367375                CPHA = CPHA + CMIN !we dont allow |c|<1m/s
    368376                ! Absolute frequency is imposed
    369                 ZO(JW, II) = CPHA * ZK(JW, II)
     377                ZO(II, JW) = CPHA * ZK(II,JW)
    370378                ! 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)
     379                ZO(II, JW) = ZO(II, JW) &
     380                     + ZK(II,JW) * COS(ZP(II,JW)) * UH(II, LAUNCH) &
     381                     + ZK(II,JW) * SIN(ZP(II,JW)) * VH(II, LAUNCH)
    374382                ! Momentum flux at launch lev
    375383                ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE
    376384                !  RIGHT IN THE SH (GWD4 after 1990)
    377                   RUW0(JW, II) = 0.
     385                  RUW0(II, JW) = 0.
    378386                 DO JJ = 1, NA
    379                     RUW0(JW, II) = RUW0(JW,II) + &
     387                    RUW0(II, JW) = RUW0(II, JW) + &
    380388         2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.)
    381389                END DO
    382                 RUW0(JW, II) = RUWFRT &
    383                           * EXP(RUW0(JW,II))/1250. &  ! 2 mpa at south pole
     390                RUW0(II,JW) = RUWFRT &
     391                          * EXP(RUW0(II,JW))/1250. &  ! 2 mpa at south pole
    384392       *((1.05+SIN(PLAT(II)*RPI/180.))/(1.01+SIN(PLAT(II)*RPI/180.))-2.05/2.01)
    385                 ! RUW0(JW, II) = RUWFRT
     393                ! RUW0(II,JW) = RUWFRT
    386394             ENDDO
    387395    end DO
    … …  
    397405
    398406       ! 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)
     407       ZOP(:, JW) = ZO(:, JW) &
     408            - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LAUNCH) &
     409            - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LAUNCH)
    402410
    403411       ! VERSION WITH FRONTAL SOURCES
    … …  
    406414
    407415       ! tanh limitation for values above CORIO (inertial instability).
    408        ! WWP(JW, :) = RUW0(JW, :) &
    409        WWP(JW, :) = RUWFRT      &
     416       ! WWP(:,JW) = RUW0(:,JW) &
     417       WWP(:,JW) = RUWFRT      &
    410418       !     * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 &
    411419       !    * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) &
    … …  
    413421       !    * (CORIO(:)*CORIO(:)) &
    414422       ! 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) &
     423       !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(:,JW)),ZOISEC)**2 &
     424       !      *ZK(:,JW)**2*DZ**2) &
    417425       ! COMPLETE FORMULA:
    418426            !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) &
    … …  
    420428       !  RESTORE DIMENSION OF A FLUX
    421429       !     *RD*TR/PR
    422        !     *1. + RUW0(JW, :)
     430       !     *1. + RUW0(:,JW)
    423431             *1.
    424432
    … …  
    429437       ! Put the stress in the right direction:
    430438
    431         RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    432         RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     439        RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     440        RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    433441
    434442    end DO
    … …  
    440448       RVW(:, LL) = 0
    441449       DO JW = 1, NW
    442           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    443           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
     450          RUW(:, LL) = RUW(:, LL) + RUWP(:,JW)
     451          RVW(:, LL) = RVW(:, LL) + RVWP(:,JW)
    444452       end DO
    445453    end DO
    … …  
    453461       ! to the next)
    454462       DO JW = 1, NW
    455           ZOM(JW, :) = ZOP(JW, :)
    456           WWM(JW, :) = WWP(JW, :)
     463          ZOM(:, JW) = ZOP(:,JW)
     464          WWM(:,JW) = WWP(:,JW)
    457465          ! Intrinsic Frequency
    458           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
    459                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
     466          ZOP(:,JW) = ZO(:, JW) - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LL + 1) &
     467               - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LL + 1)
    460468
    461469          ! No breaking (Eq.6)
    462470          ! Dissipation (Eq. 8)
    463           WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
     471          WWP(:,JW) = WWM(:,JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
    464472               + 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)))
     473               / MAX(ABS(ZOP(:,JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
     474               * ZK(:,JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
    467475
    468476          ! Critical levels (forced to zero if intrinsic frequency changes sign)
    469477          ! Saturation (Eq. 12)
    470           WWP(JW, :) = min(WWP(JW, :), MAX(0., &
    471                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
     478          WWP(:,JW) = min(WWP(:,JW), MAX(0., &
     479               SIGN(1., ZOP(:,JW) * ZOM(:, JW))) * ABS(ZOP(:,JW))**3 &
    472480          !    / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 &
    473481               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 &
    474482!              *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 &
    475483               *SATFRT**2       &
    476                / ZK(JW, :)**4)
     484               / ZK(:,JW)**4)
    477485       end DO
    478486
    … …  
    481489
    482490       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, :)
     491          RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     492          RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    485493       end DO
    486494
    … …  
    489497
    490498       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)
     499          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:,JW)
     500          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:,JW)
     501          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:,JW))/REAL(NW)
     502          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:,JW))/REAL(NW)
    495503       end DO
    496504    end DO
Note: See TracChangeset for help on using the changeset viewer.