Changeset 1283


Ignore:
Timestamp:
Jun 5, 2014, 12:34:37 PM (11 years ago)
Author:
emillour
Message:

Generic GCM:
Bug fixes:

  • hice() in physiq.F90 must be a saved array.
  • bad use of min/max on arrays in h2o_cloudrad (radii_mod.F90) which actually ended up setting all array elements to the same value.

And some cosmetic cleanup in rain.F90, vdif_kc.F and turbdiff.F90
EM

Location:
trunk/LMDZ.GENERIC
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r1216 r1283  
    10241024  over the whole planet. This should be further imroved with computation of
    10251025  means (possibly area weighed), etc.
     1026
     1027== 05/06/2014 == EM
     1028Bug fixes:
     1029- hice() in physiq.F90 must be a saved array.
     1030- bad use of min/max on arrays in h2o_cloudrad (radii_mod.F90) which actually
     1031  ended up setting all array elements to the same value.
     1032And some cosmetic cleanup in rain.F90, vdif_kc.F and turbdiff.F90
     1033

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r1252 r1283  
    348348
    349349!     included by BC for hydrology
    350       real hice(ngrid)
     350      real,allocatable,save :: hice(:)
    351351
    352352!     included by RW to test water conservation (by routine)
    … …  
    454454        ALLOCATE(cloudfrac(ngrid,nlayermx))
    455455        ALLOCATE(totcloudfrac(ngrid))
     456        ALLOCATE(hice(ngrid))
    456457        ALLOCATE(qsurf_hist(ngrid,nq))
    457458        ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
    … …  
    976977         zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
    977978         if (tracer) then
    978            pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq) 
     979           pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
    979980           dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
    980981         end if ! of if (tracer)
    … …  
    10641065         pdt(1:ngrid,1:nlayermx)    = pdt(1:ngrid,1:nlayermx) + zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx)
    10651066         zdtadj(1:ngrid,1:nlayermx) = zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
    1066  
     1067
    10671068         if(tracer) then
    10681069            pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqadj(1:ngrid,1:nlayermx,1:nq)
    … …  
    11121113              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
    11131114              zdqc)
    1114 
    11151115
    11161116         pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtc(1:ngrid,1:nlayermx)
    … …  
    12461246               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
    12471247               rainout(1:ngrid)             = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic   
    1248 
    12491248
    12501249               !-------------------------

  • trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90

    r1026 r1283  
    3838      Implicit none
    3939
    40 #include "callkeys.h"
    41 #include "dimensions.h"
    42 #include "dimphys.h"
     40      include "callkeys.h"
     41      include "dimensions.h"
     42      include "dimphys.h"
    4343
    4444      integer,intent(in) :: ngrid
    … …  
    140140      Implicit none
    141141
    142 #include "callkeys.h"
    143 #include "dimensions.h"
    144 #include "dimphys.h"
    145 #include "comcstfi.h"
     142      include "callkeys.h"
     143      include "dimensions.h"
     144      include "dimphys.h"
     145      include "comcstfi.h"
    146146
    147147      integer,intent(in) :: ngrid
    … …  
    200200      Implicit none
    201201
    202 #include "callkeys.h"
    203 #include "dimensions.h"
    204 #include "dimphys.h"
    205 #include "comcstfi.h"
     202      include "callkeys.h"
     203      include "dimensions.h"
     204      include "dimphys.h"
     205      include "comcstfi.h"
    206206
    207207      integer,intent(in) :: ngrid
    208208
    209209      real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
    210       real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (K)
     210      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (m)
    211211
    212212      real,external :: CBRT           
    213      
     213      integer :: i,k
    214214
    215215      if (radfixed) then
    … …  
    217217         reffice(1:ngrid,1:nlayermx)= rad_h2o_ice
    218218      else
    219          reffliq(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
    220          reffliq(1:ngrid,1:nlayermx)  = min(max(reffliq(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
    221          reffice(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
    222          reffice(1:ngrid,1:nlayermx)  = min(max(reffice(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
    223       end if
     219         do k=1,nlayermx
     220           do i=1,ngrid
     221             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
     222             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
     223           
     224             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
     225             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
     226           enddo
     227         enddo
     228      endif
    224229
    225230   end subroutine h2o_cloudrad
    … …  
    243248      Implicit none
    244249
    245 #include "callkeys.h"
    246 #include "dimensions.h"
    247 #include "dimphys.h"
    248 #include "comcstfi.h"
     250      include "callkeys.h"
     251      include "dimensions.h"
     252      include "dimphys.h"
     253      include "comcstfi.h"
    249254
    250255      integer,intent(in) :: ngrid,nq
    … …  
    289294      Implicit none
    290295
    291 #include "callkeys.h"
    292 #include "dimensions.h"
    293 #include "dimphys.h"
     296      include "callkeys.h"
     297      include "dimensions.h"
     298      include "dimphys.h"
    294299
    295300      integer,intent(in) :: ngrid
    … …  
    317322      Implicit none
    318323
    319 #include "callkeys.h"
    320 #include "dimensions.h"
    321 #include "dimphys.h"
     324      include "callkeys.h"
     325      include "dimensions.h"
     326      include "dimphys.h"
    322327
    323328      integer,intent(in) :: ngrid
    … …  
    346351     Implicit none
    347352
    348 #include "callkeys.h"
    349 #include "dimensions.h"
    350 #include "dimphys.h"
     353     include "callkeys.h"
     354     include "dimensions.h"
     355     include "dimphys.h"
    351356
    352357      integer,intent(in) :: ngrid

  • trunk/LMDZ.GENERIC/libf/phystd/rain.F90

    r1016 r1283  
    22
    33
    4 ! to use  'getin'
    5   use ioipsl_getincom
     4  use ioipsl_getincom, only: getin
    65  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
    76  use radii_mod, only: h2o_cloudrad
    8   USE tracer_h
     7  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    98  implicit none
    109
    … …  
    2322!==================================================================
    2423
    25 #include "dimensions.h"
    26 #include "dimphys.h"
    27 #include "comcstfi.h"
    28 #include "callkeys.h"
     24  include "dimensions.h"
     25  include "dimphys.h"
     26  include "comcstfi.h"
     27  include "callkeys.h"
    2928
    3029!     Arguments
    … …  
    149148
    150149         firstcall = .false.
    151       ENDIF
     150      ENDIF ! of IF (firstcall)
    152151
    153152!     GCM -----> subroutine variables
    … …  
    173172
    174173!     Initialise the outputs
    175       DO k = 1, nlayermx
    176       DO i = 1, ngrid
    177          d_t(i,k) = 0.0
    178          d_q(i,k) = 0.0
    179          d_ql(i,k) = 0.0
    180       ENDDO
    181       ENDDO
    182       DO i = 1, ngrid
    183          zrfl(i)  = 0.0
    184          zrfln(i) = 0.0
    185       ENDDO
     174      d_t(1:ngrid,1:nlayermx) = 0.0
     175      d_q(1:ngrid,1:nlayermx) = 0.0
     176      d_ql(1:ngrid,1:nlayermx) = 0.0
     177      zrfl(1:ngrid) = 0.0
     178      zrfln(1:ngrid) = 0.0
    186179
    187180      ! calculate saturation mixing ratio
    … …  
    203196      ENDDO
    204197
    205 
    206198      ! Vertical loop (from top to bottom)
    207199      ! We carry the rain with us and calculate that added by warm/cold precipitation
    208200      ! processes and that subtracted by evaporation at each level.
    209       DO 9999 k = nlayermx, 1, -1
     201      DO k = nlayermx, 1, -1
    210202
    211203         IF (evap_prec) THEN ! note no rneb dependence!
    … …  
    236228                     
    237229
    238                ENDIF
     230               ENDIF ! of IF (zrfl(i) .GT.0.)
    239231            ENDDO
    240          ENDIF
    241 
    242          DO i = 1, ngrid
    243             zoliq(i) = 0.0
    244          ENDDO
     232         ENDIF ! of IF (evap_prec)
     233
     234         zoliq(1:ngrid) = 0.0
    245235
    246236
    … …  
    367357         endif ! if precip_scheme=1
    368358
    369  9999 continue
     359      ENDDO ! of DO k = nlayermx, 1, -1
    370360
    371361!     Rain or snow on the ground
    … …  
    385375
    386376!     now subroutine -----> GCM variables
    387       DO k = 1, nlayermx
    388          DO i = 1, ngrid
    389            
    390             if(evap_prec)then
    391                dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
    392                d_t(i,k)          = d_t(i,k)/ptimestep
    393             else
    394                dqrain(i,k,i_vap) = 0.0
    395                d_t(i,k)          = 0.0
    396             endif
    397             dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
    398          
    399          ENDDO
    400       ENDDO
    401 
    402       RETURN
     377      if (evap_prec) then
     378        dqrain(1:ngrid,1:nlayermx,i_vap)=dqrain(1:ngrid,1:nlayermx,i_vap)/ptimestep
     379        d_t(1:ngrid,1:nlayermx)=d_t(1:ngrid,1:nlayermx)/ptimestep
     380      else
     381        dqrain(1:ngrid,1:nlayermx,i_vap)=0.0
     382        d_t(1:ngrid,1:nlayermx)=0.0
     383      endif
     384      dqrain(1:ngrid,1:nlayermx,i_ice) = d_ql(1:ngrid,1:nlayermx)/ptimestep
     385
    403386    end subroutine rain

  • trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90

    r1194 r1283  
    99      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
    1010      use radcommon_h, only : sigma, glat
    11       USE surfdat_h
    12       USE comgeomfi_h
    13       USE tracer_h
     11      use surfdat_h, only: dryness
     12      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    1413
    1514      implicit none
    … …  
    4241!     ------------
    4342
    44 #include "dimensions.h"
    45 #include "dimphys.h"
    46 #include "comcstfi.h"
    47 #include "callkeys.h"
     43      include "dimensions.h"
     44      include "dimphys.h"
     45      include "comcstfi.h"
     46      include "callkeys.h"
    4847
    4948!     arguments
    5049!     ---------
    51       INTEGER ngrid,nlay
    52       REAL ptimestep
    53       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    54       REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
    55       REAL pu(ngrid,nlay),pv(ngrid,nlay)
    56       REAL pt(ngrid,nlay),ppopsk(ngrid,nlay)
    57       REAL ptsrf(ngrid),pemis(ngrid)
    58       REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
    59       REAL pdtfi(ngrid,nlay)
    60       REAL pfluxsrf(ngrid)
    61       REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
    62       REAL pdtdif(ngrid,nlay)
    63       REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
    64       REAL pq2(ngrid,nlay+1)
     50      INTEGER,INTENT(IN) :: ngrid,nlay
     51      REAL,INTENT(IN) :: ptimestep
     52      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
     53      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
     54      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
     55      REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
     56      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
     57      REAL,INTENT(IN) :: pemis(ngrid)
     58      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
     59      REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
     60      REAL,INTENT(IN) :: pfluxsrf(ngrid)
     61      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
     62      REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
     63      REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
     64      REAL,INTENT(OUT) :: sensibFlux(ngrid)
     65      REAL,INTENT(IN) :: pcapcal(ngrid)
     66      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
    6567     
    66       integer rnat(ngrid)     
     68      integer,intent(in) :: rnat(ngrid)     
     69      LOGICAL,INTENT(IN) :: lastcall ! not used
    6770
    6871!     Arguments added for condensation
    69       logical lecrit
    70       REAL pz0
     72      logical,intent(in) :: lecrit ! not used.
     73      REAL,INTENT(IN) :: pz0
    7174
    7275!     Tracers
    7376!     --------
    74       integer nq
    75       real pqsurf(ngrid,nq)
    76       real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
    77       real pdqdif(ngrid,nlay,nq)
    78       real pdqsdif(ngrid,nq)
     77      integer,intent(in) :: nq
     78      real,intent(in) :: pqsurf(ngrid,nq)
     79      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
     80      real,intent(out) :: pdqdif(ngrid,nlay,nq)
     81      real,intent(out) :: pdqsdif(ngrid,nq)
    7982     
    8083!     local
    … …  
    102105      REAL zx_alf1(ngrid),zx_alf2(ngrid)
    103106
    104       LOGICAL firstcall
    105       SAVE firstcall
     107      LOGICAL,SAVE :: firstcall=.true.
    106108     
    107       LOGICAL lastcall
    108109!     Tracers
    109110!     -------
    … …  
    115116      REAL rho(ngrid)         ! near-surface air density
    116117      REAL qsat(ngrid),psat(ngrid)
    117       DATA firstcall/.true./
    118118      REAL kmixmin
    119119
    … …  
    242242
    243243      call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
    244 
     244     
    245245!     Adding eddy mixing to mimic 3D general circulation in 1D
    246246!     R. Wordsworth & F. Forget (2010)
    … …  
    731731!      endif
    732732
    733 
    734       return
    735733      end

  • trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F

    r787 r1283  
    2727c
    2828c.......................................................................
    29       INTEGER ngrid
    30       REAL dt,g
    31       REAL zlev(ngrid,nlayermx+1)
    32       REAL zlay(ngrid,nlayermx)
    33       REAL u(ngrid,nlayermx)
    34       REAL v(ngrid,nlayermx)
    35       REAL teta(ngrid,nlayermx)
    36       REAL cd(ngrid)
    37       REAL q2(ngrid,nlayermx+1)
    38       REAL km(ngrid,nlayermx+1)
    39       REAL kn(ngrid,nlayermx+1)
     29      INTEGER,INTENT(IN) :: ngrid
     30      REAL,INTENT(IN) :: dt,g
     31      REAL,INTENT(IN) :: zlev(ngrid,nlayermx+1)
     32      REAL,INTENT(IN) :: zlay(ngrid,nlayermx)
     33      REAL,INTENT(IN) :: u(ngrid,nlayermx)
     34      REAL,INTENT(IN) :: v(ngrid,nlayermx)
     35      REAL,INTENT(IN) :: teta(ngrid,nlayermx)
     36      REAL,INTENT(IN) :: cd(ngrid)
     37      REAL,INTENT(INOUT) :: q2(ngrid,nlayermx+1)
     38      REAL,INTENT(OUT) :: km(ngrid,nlayermx+1)
     39      REAL,INTENT(OUT) :: kn(ngrid,nlayermx+1)
    4040c.......................................................................
    4141c
    … …  
    5050c
    5151c.......................................................................
    52       INTEGER nlay,nlev
     52      INTEGER,PARAMETER :: nlay=nlayermx
     53      INTEGER,PARAMETER :: nlev=nlayermx+1
    5354      REAL unsdz(ngrid,nlayermx)
    5455      REAL unsdzdec(ngrid,nlayermx+1)
    … …  
    114115c.......................................................................
    115116      REAL gn
    116       REAL gnmin
    117       REAL gnmax
     117      REAL,PARAMETER :: gnmin=-10.E+0
     118      REAL,PARAMETER :: gnmax=0.0233E+0
    118119      LOGICAL gninf
    119120      LOGICAL gnsup
    … …  
    134135c
    135136c.......................................................................
    136       REAL kappa
    137       REAL long0
    138       REAL a1,a2,b1,b2,c1
    139       REAL cn1,cn2
    140       REAL cm1,cm2,cm3,cm4
     137      REAL,PARAMETER :: kappa=0.4E+0
     138      REAL,PARAMETER :: long0=160.E+0
     139      REAL,PARAMETER :: a1=0.92E+0,a2=0.74E+0
     140      REAL,PARAMETER :: b1=16.6E+0,b2=10.1E+0,c1=0.08E+0
     141      REAL,PARAMETER :: cn1=a2*(1.E+0 -6.E+0 *a1/b1)
     142      REAL,PARAMETER :: cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
     143      REAL,PARAMETER :: cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
     144      REAL,PARAMETER :: cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*
     145     &          (1.E+0 -6.E+0 *a1/b1)-3.E+0 *c1*(b2+6.E+0 *a1)))
     146      REAL,PARAMETER :: cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
     147      REAL,PARAMETER :: cm4=-9.E+0 *a1*a2
    141148c.......................................................................
    142149c
    … …  
    157164c
    158165c.......................................................................
    159       REAL q2min
    160       REAL q2max
     166      REAL,PARAMETER :: q2min=1.E-3
     167      REAL,PARAMETER :: q2max=1.E+2
    161168c.......................................................................
    162169c knmin : borne inferieure de kn
    163170c kmmin : borne inferieure de km
    164171c.......................................................................
    165       REAL knmin
    166       REAL kmmin
     172      REAL,PARAMETER :: knmin=1.E-5
     173      REAL,PARAMETER :: kmmin=1.E-5
    167174c.......................................................................
    168175      INTEGER ilay,ilev,igrid
    169176      REAL tmp1,tmp2
    170177c.......................................................................
    171       PARAMETER (kappa=0.4E+0)
    172       PARAMETER (long0=160.E+0)
    173       PARAMETER (gnmin=-10.E+0)
    174       PARAMETER (gnmax=0.0233E+0)
    175       PARAMETER (a1=0.92E+0)
    176       PARAMETER (a2=0.74E+0)
    177       PARAMETER (b1=16.6E+0)
    178       PARAMETER (b2=10.1E+0)
    179       PARAMETER (c1=0.08E+0)
    180       PARAMETER (knmin=1.E-5)
    181       PARAMETER (kmmin=1.E-5)
    182       PARAMETER (q2min=1.E-3)
    183       PARAMETER (q2max=1.E+2)
    184       PARAMETER (nlay=nlayermx)
    185       PARAMETER (nlev=nlayermx+1)
    186 c
    187       PARAMETER (
    188      &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
    189      &          )
    190       PARAMETER (
    191      &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
    192      &          )
    193       PARAMETER (
    194      &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
    195      &          )
    196       PARAMETER (
    197      &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
    198      &          -3.E+0 *c1*(b2+6.E+0 *a1)))
    199      &          )
    200       PARAMETER (
    201      &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
    202      &          )
    203       PARAMETER (
    204      &  cm4=-9.E+0 *a1*a2
    205      &          )
     178c
     179
     180! initialization of local variables:
     181      long(:,:)=0.
     182      n2(:,:)=0.
     183      sn(:,:)=0.
     184      snq2(:,:)=0.
     185      sm(:,:)=0.
     186      smq2(:,:)=0.
     187
    206188c.......................................................................
    207189c  traitment des valeur de q2 en entree
    … …  
    209191c
    210192      DO ilev=1,nlev
    211                                                  DO igrid=1,ngrid
     193       DO igrid=1,ngrid
    212194        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
    213195        q(igrid,ilev)=sqrt(q2(igrid,ilev))
    214                                                  ENDDO
    215       ENDDO
    216 c
    217                                                  DO igrid=1,ngrid
    218       tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
    219       q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
    220       q2(igrid,1)=amax1(q2(igrid,1),q2min)
    221       q(igrid,1)=sqrt(q2(igrid,1))
    222                                                  ENDDO
     196       ENDDO
     197      ENDDO
     198c
     199      DO igrid=1,ngrid
     200        tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
     201        q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
     202        q2(igrid,1)=amax1(q2(igrid,1),q2min)
     203        q(igrid,1)=sqrt(q2(igrid,1))
     204      ENDDO
    223205c
    224206c.......................................................................
    … …  
    237219c
    238220      DO ilay=1,nlay
    239                                                  DO igrid=1,ngrid
     221       DO igrid=1,ngrid
    240222        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
    241                                                  ENDDO
    242       ENDDO
    243                                                  DO igrid=1,ngrid
    244       unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
    245                                                  ENDDO
     223       ENDDO
     224      ENDDO
     225
     226      DO igrid=1,ngrid
     227        unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
     228      ENDDO
     229
    246230      DO ilay=2,nlay
    247                                                  DO igrid=1,ngrid
    248         unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
    249                                                  ENDDO
    250       ENDDO
    251                                                  DO igrid=1,ngrid
    252       unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
    253                                                  ENDDO
     231        DO igrid=1,ngrid
     232          unsdzdec(igrid,ilay)=1.E+0/
     233     &                          (zlay(igrid,ilay)-zlay(igrid,ilay-1))
     234        ENDDO
     235      ENDDO
     236     
     237      DO igrid=1,ngrid
     238        unsdzdec(igrid,nlay+1)=1.E+0/
     239     &                          (zlev(igrid,nlay+1)-zlay(igrid,nlay))
     240      ENDDO
    254241c
    255242c.......................................................................
    … …  
    257244c.......................................................................
    258245c
    259                                                  DO igrid=1,ngrid
    260       m2(igrid,1)=(unsdzdec(igrid,1)
     246      DO igrid=1,ngrid
     247        m2(igrid,1)=(unsdzdec(igrid,1)
    261248     &                   *u(igrid,1))**2
    262249     &                 +(unsdzdec(igrid,1)
    263250     &                   *v(igrid,1))**2
    264       m(igrid,1)=sqrt(m2(igrid,1))
    265       mpre(igrid,1)=m(igrid,1)
    266                                                  ENDDO
     251        m(igrid,1)=sqrt(m2(igrid,1))
     252        mpre(igrid,1)=m(igrid,1)
     253      ENDDO
    267254c
    268255c-----------------------------------------------------------------------
    269256      DO ilev=2,nlev-1
    270                                                  DO igrid=1,ngrid
     257       DO igrid=1,ngrid
    271258c-----------------------------------------------------------------------
    272259c
    … …  
    295282c
    296283c-----------------------------------------------------------------------
    297                                                  ENDDO
    298       ENDDO
    299 c-----------------------------------------------------------------------
    300 c
    301                                                  DO igrid=1,ngrid
    302       m2(igrid,nlev)=m2(igrid,nlev-1)
    303       m(igrid,nlev)=m(igrid,nlev-1)
    304       mpre(igrid,nlev)=m(igrid,nlev)
    305                                                  ENDDO
     284       ENDDO
     285      ENDDO
     286c-----------------------------------------------------------------------
     287c
     288      DO igrid=1,ngrid
     289        m2(igrid,nlev)=m2(igrid,nlev-1)
     290        m(igrid,nlev)=m(igrid,nlev-1)
     291        mpre(igrid,nlev)=m(igrid,nlev)
     292      ENDDO
     293     
    306294c
    307295c.......................................................................
    … …  
    311299c-----------------------------------------------------------------------
    312300      DO ilev=2,nlev-1
    313                                                  DO igrid=1,ngrid
     301       DO igrid=1,ngrid
    314302c-----------------------------------------------------------------------
    315303c
    … …  
    375363c
    376364c-----------------------------------------------------------------------
    377                                                   ENDDO
    378       ENDDO
     365       ENDDO ! of DO igrid=1,ngrid
     366      ENDDO ! of DO ilev=2,nlev-1
    379367c-----------------------------------------------------------------------
    380368c
    … …  
    383371c.......................................................................
    384372c
    385                                                  DO igrid=1,ngrid
    386       kn(igrid,1)=knmin
    387       km(igrid,1)=kmmin
    388       kmpre(igrid,1)=km(igrid,1)
    389                                                  ENDDO
     373      DO igrid=1,ngrid
     374        kn(igrid,1)=knmin
     375        km(igrid,1)=kmmin
     376        kmpre(igrid,1)=km(igrid,1)
     377      ENDDO
    390378c
    391379c-----------------------------------------------------------------------
    392380      DO ilev=2,nlev-1
    393                                                  DO igrid=1,ngrid
    394 c-----------------------------------------------------------------------
    395 c
     381       DO igrid=1,ngrid
    396382        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
    397383     &                                         *sn(igrid,ilev)
    … …  
    399385     &                                         *sm(igrid,ilev)
    400386        kmpre(igrid,ilev)=km(igrid,ilev)
    401 c
    402 c-----------------------------------------------------------------------
    403                                                  ENDDO
    404       ENDDO
    405 c-----------------------------------------------------------------------
    406 c
    407                                                  DO igrid=1,ngrid
    408       kn(igrid,nlev)=kn(igrid,nlev-1)
    409       km(igrid,nlev)=km(igrid,nlev-1)
    410       kmpre(igrid,nlev)=km(igrid,nlev)
    411                                                  ENDDO
     387       ENDDO
     388      ENDDO
     389c-----------------------------------------------------------------------
     390c
     391      DO igrid=1,ngrid
     392        kn(igrid,nlev)=kn(igrid,nlev-1)
     393        km(igrid,nlev)=km(igrid,nlev-1)
     394        kmpre(igrid,nlev)=km(igrid,nlev)
     395      ENDDO
    412396c
    413397c.......................................................................
    … …  
    539523c
    540524c
    541                                                  DO igrid=1,ngrid
    542       kn(igrid,1)=knmin
    543       km(igrid,1)=kmmin
    544       q2(igrid,nlev)=q2(igrid,nlev-1)
    545       q(igrid,nlev)=q(igrid,nlev-1)
    546       kn(igrid,nlev)=kn(igrid,nlev-1)
    547       km(igrid,nlev)=km(igrid,nlev-1)
    548                                                  ENDDO
     525      DO igrid=1,ngrid
     526        kn(igrid,1)=knmin
     527        km(igrid,1)=kmmin
     528        q2(igrid,nlev)=q2(igrid,nlev-1)
     529        q(igrid,nlev)=q(igrid,nlev-1)
     530        kn(igrid,nlev)=kn(igrid,nlev-1)
     531        km(igrid,nlev)=km(igrid,nlev-1)
     532      ENDDO
    549533
    550 c
    551 c.......................................................................
    552 c
    553       RETURN
    554534      END
Note: See TracChangeset for help on using the changeset viewer.