Changeset 2729


Ignore:
Timestamp:
Dec 13, 2016, 5:13:39 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2719:2727 into testing branch

Location:
LMDZ5/branches/testing
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/DefLists/field_def_lmdz.xml

    r2720 r2729  
    880880      <field id="c06_clmodis" long_name="MODIS Cloud Area Fraction" unit="1" axis_ref="pressure2" />
    881881      <field id="c07_clmodis" long_name="MODIS Cloud Area Fraction" unit="1" axis_ref="pressure2" />
     882      <field id="c01_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     883      <field id="c02_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     884      <field id="c03_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     885      <field id="c04_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     886      <field id="c05_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     887      <field id="c06_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     888      <field id="c07_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" />
     889      <field id="c01_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     890      <field id="c02_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     891      <field id="c03_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     892      <field id="c04_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     893      <field id="c05_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     894      <field id="c06_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
     895      <field id="c07_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" />
    882896    </field_group>
    883897
  • LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90

    r2542 r2729  
    5656       ENDIF
    5757
    58        IF (size(time).NE.year_len) THEN
     58!--test if time is different from year_len but allow a mismatch of 1 day
     59       IF (size(time).NE.year_len.AND.size(time).NE.year_len+1) THEN
    5960         PRINT *,'read_rsun_rrtm time <> year_len = ', size(time), year_len
    6061         CALL abort_physic('read_rsun_rrtm','time dim should be the number of days in year',1)
     62       ENDIF
     63!--warning only if forcing file has 366 days but year_len has only 365
     64       IF (size(time).EQ.year_len+1) THEN
     65         PRINT *,'Warning read_rsun_rrtm uses a leap year rsun for a noleap year'
    6166       ENDIF
    6267
     
    9196
    9297!--copy
    93       RSUN(:)=SSI_FRAC(:,day_cur)
     98      RSUN(1:NSW)=SSI_FRAC(:,day_cur)
    9499      solaire=TSI(day_cur)
    95100
  • LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90

    r2720 r2729  
    8181
    8282!--only root reads the data
    83     IF (is_mpi_root.AND.is_omp_root) THEN
     83      IF (is_mpi_root.AND.is_omp_root) THEN
    8484
    8585!--check mth_cur
    86     IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
    87       print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
    88     ENDIF
     86        IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
     87          print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
     88        ENDIF
    8989
    9090!--initialize n_lon as input data is 2D (lat-alt) only
    91     n_lon = nbp_lon
     91        n_lon = nbp_lon
    9292
    9393!--Starts with SW optical properties
    9494
    95     CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
    96 
    97     CALL nf95_inq_varid(ncid_in, "LEV", varid)
    98     CALL nf95_gw_var(ncid_in, varid, lev)
    99     n_lev = size(lev)
    100     IF (n_lev.NE.klev) THEN
    101        abort_message='Le nombre de niveaux n est pas egal a klev'
    102        CALL abort_physic(modname,abort_message,1)
    103     ENDIF
    104 
    105     CALL nf95_inq_varid(ncid_in, "LAT", varid)
    106     CALL nf95_gw_var(ncid_in, varid, latitude)
    107     n_lat = size(latitude)
    108     IF (n_lat.NE.nbp_lat) THEN
    109        print *, 'latitude=', n_lat, nbp_lat
    110        abort_message='Le nombre de lat n est pas egal a nbp_lat'
    111        CALL abort_physic(modname,abort_message,1)
    112     ENDIF
    113 
    114     CALL nf95_inq_varid(ncid_in, "TIME", varid)
    115     CALL nf95_gw_var(ncid_in, varid, time)
    116     n_month = size(time)
    117     IF (n_month.NE.12) THEN
    118        abort_message='Le nombre de month n est pas egal a 12'
    119        CALL abort_physic(modname,abort_message,1)
    120     ENDIF
    121 
    122     CALL nf95_inq_varid(ncid_in, "WAV", varid)
    123     CALL nf95_gw_var(ncid_in, varid, wav)
    124     n_wav = size(wav)
    125     print *, 'WAV aerosol strato=', n_wav, wav
    126     IF (n_wav.NE.NSW) THEN
    127        abort_message='Le nombre de wav n est pas egal a NSW'
    128        CALL abort_physic(modname,abort_message,1)
    129     ENDIF
    130 
    131     ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
    132     ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
    133     ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
    134 
    135     ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    136     ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    137     ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    138 
    139     ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    140     ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    141     ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     95        CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
     96
     97        CALL nf95_inq_varid(ncid_in, "LEV", varid)
     98        CALL nf95_gw_var(ncid_in, varid, lev)
     99        n_lev = size(lev)
     100        IF (n_lev.NE.klev) THEN
     101           abort_message='Le nombre de niveaux n est pas egal a klev'
     102           CALL abort_physic(modname,abort_message,1)
     103        ENDIF
     104
     105        CALL nf95_inq_varid(ncid_in, "LAT", varid)
     106        CALL nf95_gw_var(ncid_in, varid, latitude)
     107        n_lat = size(latitude)
     108        IF (n_lat.NE.nbp_lat) THEN
     109           print *, 'latitude=', n_lat, nbp_lat
     110           abort_message='Le nombre de lat n est pas egal a nbp_lat'
     111           CALL abort_physic(modname,abort_message,1)
     112        ENDIF
     113
     114        CALL nf95_inq_varid(ncid_in, "TIME", varid)
     115        CALL nf95_gw_var(ncid_in, varid, time)
     116        n_month = size(time)
     117        IF (n_month.NE.12) THEN
     118           abort_message='Le nombre de month n est pas egal a 12'
     119           CALL abort_physic(modname,abort_message,1)
     120        ENDIF
     121
     122        CALL nf95_inq_varid(ncid_in, "WAV", varid)
     123        CALL nf95_gw_var(ncid_in, varid, wav)
     124        n_wav = size(wav)
     125        print *, 'WAV aerosol strato=', n_wav, wav
     126        IF (n_wav.NE.NSW) THEN
     127           abort_message='Le nombre de wav n est pas egal a NSW'
     128           CALL abort_physic(modname,abort_message,1)
     129        ENDIF
     130
     131        ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
     132        ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
     133        ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
     134
     135        ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     136        ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     137        ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     138
     139        ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     140        ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     141        ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    142142
    143143!--reading stratospheric aerosol tau per layer
    144     CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
    145     ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
    146     print *,'code erreur readaerosolstrato=', ncerr, varid
     144        CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
     145        ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
     146        print *,'code erreur readaerosolstrato=', ncerr, varid
    147147
    148148!--reading stratospheric aerosol omega per layer
    149     CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
    150     ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
    151     print *,'code erreur readaerosolstrato=', ncerr, varid
     149        CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
     150        ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
     151        print *,'code erreur readaerosolstrato=', ncerr, varid
    152152
    153153!--reading stratospheric aerosol g per layer
    154     CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
    155     ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
    156     print *,'code erreur readaerosolstrato sw=', ncerr, varid
    157 
    158     CALL nf95_close(ncid_in)
     154        CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
     155        ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
     156        print *,'code erreur readaerosolstrato sw=', ncerr, varid
     157
     158        CALL nf95_close(ncid_in)
    159159
    160160!--select the correct month
    161161!--and copy into 1st longitude
    162     tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
    163     pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
    164     cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
     162        tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
     163        pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
     164        cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
    165165
    166166!--copy longitudes
    167     DO i=2, n_lon
    168      tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
    169      pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
    170       cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
    171     ENDDO
     167        DO i=2, n_lon
     168         tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
     169         pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
     170         cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
     171        ENDDO
    172172
    173173!---reduce to a klon_glo grid
    174     DO band=1, NSW
    175       CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
    176       CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
    177       CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
    178     ENDDO
     174        DO band=1, NSW
     175          CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
     176          CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
     177          CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
     178        ENDDO
    179179
    180180!--Now LW optical properties
    181181!
    182     CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
    183 
    184     CALL nf95_inq_varid(ncid_in, "LEV", varid)
    185     CALL nf95_gw_var(ncid_in, varid, lev)
    186     n_lev = size(lev)
    187     IF (n_lev.NE.klev) THEN
    188        abort_message='Le nombre de niveaux n est pas egal a klev'
    189        CALL abort_physic(modname,abort_message,1)
    190     ENDIF
    191 
    192     CALL nf95_inq_varid(ncid_in, "LAT", varid)
    193     CALL nf95_gw_var(ncid_in, varid, latitude)
    194     n_lat = size(latitude)
    195     IF (n_lat.NE.nbp_lat) THEN
    196        abort_message='Le nombre de lat n est pas egal a nbp_lat'
    197        CALL abort_physic(modname,abort_message,1)
    198     ENDIF
    199 
    200     CALL nf95_inq_varid(ncid_in, "TIME", varid)
    201     CALL nf95_gw_var(ncid_in, varid, time)
    202     n_month = size(time)
    203     IF (n_month.NE.12) THEN
    204        abort_message='Le nombre de month n est pas egal a 12'
    205        CALL abort_physic(modname,abort_message,1)
    206     ENDIF
    207 
    208     CALL nf95_inq_varid(ncid_in, "WAV", varid)
    209     CALL nf95_gw_var(ncid_in, varid, wav)
    210     n_wav = size(wav)
    211     print *, 'WAV aerosol strato=', n_wav, wav
    212     IF (n_wav.NE.NLW) THEN
    213        abort_message='Le nombre de wav n est pas egal a NLW'
    214        CALL abort_physic(modname,abort_message,1)
    215     ENDIF
    216 
    217     ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
    218     ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    219     ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     182        CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
     183
     184        CALL nf95_inq_varid(ncid_in, "LEV", varid)
     185        CALL nf95_gw_var(ncid_in, varid, lev)
     186        n_lev = size(lev)
     187        IF (n_lev.NE.klev) THEN
     188           abort_message='Le nombre de niveaux n est pas egal a klev'
     189           CALL abort_physic(modname,abort_message,1)
     190        ENDIF
     191
     192        CALL nf95_inq_varid(ncid_in, "LAT", varid)
     193        CALL nf95_gw_var(ncid_in, varid, latitude)
     194        n_lat = size(latitude)
     195        IF (n_lat.NE.nbp_lat) THEN
     196           abort_message='Le nombre de lat n est pas egal a nbp_lat'
     197           CALL abort_physic(modname,abort_message,1)
     198        ENDIF
     199
     200        CALL nf95_inq_varid(ncid_in, "TIME", varid)
     201        CALL nf95_gw_var(ncid_in, varid, time)
     202        n_month = size(time)
     203        IF (n_month.NE.12) THEN
     204           abort_message='Le nombre de month n est pas egal a 12'
     205           CALL abort_physic(modname,abort_message,1)
     206        ENDIF
     207
     208        CALL nf95_inq_varid(ncid_in, "WAV", varid)
     209        CALL nf95_gw_var(ncid_in, varid, wav)
     210        n_wav = size(wav)
     211        print *, 'WAV aerosol strato=', n_wav, wav
     212        IF (n_wav.NE.NLW) THEN
     213           abort_message='Le nombre de wav n est pas egal a NLW'
     214           CALL abort_physic(modname,abort_message,1)
     215        ENDIF
     216
     217        ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
     218        ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     219        ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    220220
    221221!--reading stratospheric aerosol lw tau per layer
    222     CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
    223     ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
    224     print *,'code erreur readaerosolstrato lw=', ncerr, varid
    225 
    226     CALL nf95_close(ncid_in)
     222        CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
     223        ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
     224        print *,'code erreur readaerosolstrato lw=', ncerr, varid
     225
     226        CALL nf95_close(ncid_in)
    227227
    228228!--select the correct month
    229229!--and copy into 1st longitude
    230     taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
     230        taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
    231231!--copy longitudes
    232     DO i=2, n_lon
    233       taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
    234     ENDDO
     232        DO i=2, n_lon
     233          taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
     234        ENDDO
    235235
    236236!---reduce to a klon_glo grid
    237     DO band=1, NLW
    238       CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
    239     ENDDO
    240 
    241     ENDIF !--is_mpi_root and is_omp_root
     237        DO band=1, NLW
     238          CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
     239        ENDDO
     240
     241      ELSE !--proc other than mpi_root and omp_root
     242           !--dummy allocation needed for debug mode
     243
     244        ALLOCATE(tauaerstrat_mois_glo(1,1,1))
     245        ALLOCATE(pizaerstrat_mois_glo(1,1,1))
     246        ALLOCATE(cgaerstrat_mois_glo(1,1,1))
     247        ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
     248
     249      ENDIF !--is_mpi_root and is_omp_root
    242250
    243251!$OMP BARRIER
    244252
    245253!--keep memory of previous month
    246     mth_pre=mth_cur
     254      mth_pre=mth_cur
    247255
    248256!--scatter on all proc
    249     CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
    250     CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
    251     CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
    252     CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
    253 
    254     IF (is_mpi_root.AND.is_omp_root) THEN
    255 !
    256     DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
    257     DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
    258     DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
    259 
    260     DEALLOCATE(taulwaerstrat,taulwaerstrat_mois,taulwaerstrat_mois_glo)
    261 !
    262     ENDIF !--is_mpi_root and is_omp_root
     257      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
     258      CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
     259      CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
     260      CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
     261
     262      IF (is_mpi_root.AND.is_omp_root) THEN
     263!
     264        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
     265        DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
     266        DEALLOCATE(taulwaerstrat,taulwaerstrat_mois)
     267!
     268      ENDIF !--is_mpi_root and is_omp_root
     269
     270      DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
     271      DEALLOCATE(taulwaerstrat_mois_glo)
    263272
    264273!$OMP BARRIER
     
    282291!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
    283292    DO band=1, NSW
    284     WHERE (stratomask.GT.0.999999)
     293      WHERE (stratomask.GT.0.999999)
    285294!--anthropogenic aerosols bands 1 to NSW
    286     cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
    287                                      cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
    288                                 MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
    289                                      piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
    290     piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
    291                                      piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
    292                                 MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
    293     tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
     295        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
     296                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
     297                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
     298                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
     299        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
     300                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
     301                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
     302        tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
    294303!--natural aerosols bands 1 to NSW
    295     cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
    296                                      cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
    297                                 MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
    298                                      piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
    299     piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
    300                                      piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
    301                                 MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
    302     tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
     304        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
     305                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
     306                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
     307                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
     308        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
     309                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
     310                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
     311        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
    303312!--no stratospheric aerosol in index 1 for these tests
    304 !    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
    305 !    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
    306 !    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
     313!        cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
     314!        piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
     315!        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
    307316    ENDWHERE
    308317    ENDDO
     
    322331
    323332    DO band=1, NLW
    324     WHERE (stratomask.GT.0.999999)
    325     tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
    326     tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
     333      WHERE (stratomask.GT.0.999999)
     334        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
     335        tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
    327336!--no stratospheric aerosols in index 1 for these tests
    328337!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band)
    329     ENDWHERE
     338      ENDWHERE
    330339    ENDDO
    331340
  • LMDZ5/branches/testing/libf/phylmd/yamada4.F90

    r2594 r2729  
    22
    33SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
    4     cd, q2, km, kn, kq, ustar, iflag_pbl)
     4    cd, tke, km, kn, kq, ustar, iflag_pbl)
    55
    66  USE dimphy
     
    5656  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
    5757  !                          -> the model can run with longer time-steps.
     58  !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
     59  !               On met tke (=q2/2) en entrée plutôt que q2
     60  !               On corrige l'update de la tke
    5861  !
    5962  ! Inpout/Output :
    6063  !==============
    61   ! q2 : $q^2$ au bas de chaque couche
     64  ! tke : tke au bas de chaque couche
    6265  ! (en entree : la valeur au debut du pas de temps)
    6366  ! (en sortie : la valeur a la fin du pas de temps)
     
    9093  REAL teta(klon, klev)
    9194  REAL cd(klon)
    92   REAL q2(klon, klev+1)
     95  REAL tke(klon, klev+1)
    9396  REAL unsdz(klon, klev)
    9497  REAL unsdzdec(klon, klev+1)
     
    104107  INCLUDE "clesphys.h"
    105108
    106 
     109  REAL q2(klon, klev+1)
    107110  REAL kmpre(klon, klev+1), tmp2, qpre
    108111  REAL mpre(klon, klev+1)
     
    127130  REAL zz(klon, klev+1)
    128131  INTEGER iter
    129   REAL ric, rifc, b1, kap
    130   SAVE ric, rifc, b1, kap
    131   DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
    132   !$OMP THREADPRIVATE(ric,rifc,b1,kap)
    133   REAL seuilsm, seuilalpha
     132  REAL,SAVE :: ric0,ric,rifc, b1, kap
     133  !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
     134  DATA b1, kap/16.6, 0.4/
     135  REAL,SAVE :: seuilsm, seuilalpha
     136  !$OMP THREADPRIVATE(seuilsm, seuilalpha)
    134137  REAL,SAVE :: lmixmin
    135138  !$OMP THREADPRIVATE(lmixmin)
     139  LOGICAL, SAVE :: new_yamada4
     140  !$OMP THREADPRIVATE(new_yamada4)
     141  REAL, SAVE :: yun,ydeux
     142  !$OMP THREADPRIVATE(yun,ydeux)
    136143  REAL frif, falpha, fsm
    137144  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
     
    141148
    142149
     150
    143151  ! Fonctions utilis??es
    144152  !--------------------
     
    150158
    151159  IF (firstcall) THEN
     160! Seuil dans le code de turbulence
     161    new_yamada4=.false.
     162    CALL getin_p('new_yamada4',new_yamada4)
     163! Pour garantir la continuite dans la mise au point de CMIP6.
     164! Eliminer l'option new_yamada4 des le printemps 2016
     165    IF (new_yamada4) THEN
     166       ric=0.143 ! qui donne des valeurs proches des seuils proposes
     167                 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
     168                 !                    sm=1.1213 (1.12  dans Y83)
     169       CALL getin_p('yamada4_ric',ric)
     170       ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
     171       ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif
     172       rifc=frif(ric)
     173       seuilsm=fsm(frif(ric))
     174       seuilalpha=falpha(frif(ric))
     175       yun=1.
     176       ydeux=2.
     177!      yun=2.
     178!      ydeux=1.
     179    ELSE
     180  !!  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
     181       ric=0.195
     182       rifc=0.191
     183       seuilalpha=1.12
     184       seuilsm=0.085
     185       yun=2.
     186       ydeux=1.
     187    ENDIF
     188    PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
    152189    firstcall = .FALSE.
    153190    lmixmin=1.
     
    169206  END IF
    170207
    171 ! Seuil dans le code de turbulence
    172  seuilalpha=1.12
    173  seuilsm=0.085
    174208
    175209  nlay = klev
     
    231265        rif(ig, k) = rifc
    232266      END IF
    233 
     267if (new_yamada4) then
     268        alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
     269        sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
     270else
    234271      IF (rif(ig,k)<0.16) THEN
    235272        alpha(ig, k) = falpha(rif(ig,k))
     
    239276        sm(ig, k) = seuilsm
    240277      END IF
     278
     279end if
    241280      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
    242281    END DO
     
    244283
    245284
     285
     286
     287
     288  !=======================================================================
     289  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
     290  !=======================================================================
     291
     292  ! On commence par calculé q2 à partir de la tke
     293
     294  IF (new_yamada4) THEN
     295      DO k=1,klev+1
     296         q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
     297      ENDDO
     298  ELSE
     299      DO k=1,klev+1
     300         q2(1:ngrid,k)=tke(1:ngrid,k)
     301      ENDDO
     302  ENDIF
    246303
    247304! ====================================================================
     
    252309  CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
    253310
    254 
    255 
    256   !=======================================================================
    257   !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
    258   !=======================================================================
    259311
    260312  !--------------
     
    350402    DO k = 2, klev - 1
    351403      km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
    352       q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
     404      q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
     405!     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
    353406      q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
    354       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
     407       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
     408!     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
    355409      q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
    356410    END DO
     
    517571      lyam(1:ngrid, 2:klev)
    518572  END IF
     573
     574
     575
     576!============================================================================
     577! Mise à jour de la tke
     578!============================================================================
     579
     580  IF (new_yamada4) THEN
     581     DO k=1,klev+1
     582        tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
     583     ENDDO
     584  ELSE
     585     DO k=1,klev+1
     586        tke(1:ngrid,k)=q2(1:ngrid,k)
     587     ENDDO
     588  ENDIF
     589
    519590
    520591!============================================================================
Note: See TracChangeset for help on using the changeset viewer.