Commit 1ca781e42b083260b5dbe320a1f74c923dd38416

Authored by dumoda01
1 parent a6b0c648
Exists in master and in 1 other branch snow

Il est maintenant possible de specifier des valeurs de concentration de glace po…

…ur attenuer la lumiere incidente a la surface de l'ocean a l'aide d'un fichier ice.dat. Les options sont documentees dans le fichier nml/airsea.nml qui a 3 parametres supplementaires. Il se peut donc que vous deviez absolument tenir compte de ces changements dans votre airsea.nml afin d'eviter des problemes. Les modifications sont toutes contenues dans airsea.F90 et airsea.nml.
Showing 2 changed files with 129 additions and 6 deletions   Show diff stats
nml/airsea.nml
... ... @@ -21,6 +21,19 @@
21 21 ! or dew point temperature in C (depending on wet_mode)
22 22 ! - cloud cover in 1/10
23 23 !
  24 +! ice_method -> method for taking sea ice account in light attenuation
  25 +! 0: No sea ice (ice=0)
  26 +! 1: Constant ice concentration
  27 +! 2: ice are read from ice_file
  28 +!
  29 +! ice_file -> file with sea ice data (for calc_fluxes=.true.)
  30 +! used to attenuation the shortwave radiation at sea surface
  31 +! - ice concentration in %
  32 +!
  33 +! iceatt -> Attenuation coefficient through ice
  34 +! (0.96 means that 4% of light will pass
  35 + through a complete ice cover)
  36 +!
24 37 ! wet_mode -> decides what is given in 7. column in meteo_file
25 38 ! 1: relative humidity
26 39 ! 2: wet bulb temperature
... ... @@ -79,6 +92,9 @@
79 92 &airsea
80 93 calc_fluxes= .false.
81 94 meteo_file= 'meteo.dat'
  95 + ice_method= 0
  96 + ice_file= 'ice.dat'
  97 + iceatt= 0.90
82 98 wet_mode= 1
83 99 heat_method= 2
84 100 const_swr= 100.0
... ...
src/airsea/airsea.F90
... ... @@ -55,7 +55,8 @@
55 55  
56 56 ! sea surface temperature (degC) and
57 57 ! sea surface salinity (psu)
58   - REALTYPE, public :: sst,sss
  58 +!DD sea ice concentration
  59 + REALTYPE, public :: sst,sss,ice !DD
59 60  
60 61 ! integrated short-wave radiation,
61 62 ! surface heat flux (J/m^2)
... ... @@ -72,6 +73,7 @@
72 73 integer, parameter :: p_e_unit=23
73 74 integer, parameter :: sst_unit=24
74 75 integer, parameter :: sss_unit=25
  76 + integer, parameter :: ice_unit=26
75 77  
76 78 REALTYPE, parameter :: cpa=1008.
77 79 REALTYPE, parameter :: cp=3985.
... ... @@ -148,6 +150,7 @@
148 150 integer :: p_e_method
149 151 integer :: sst_method
150 152 integer :: sss_method
  153 + integer :: ice_method !DD
151 154 integer :: wet_mode
152 155  
153 156 character(len=PATH_MAX) :: meteo_file
... ... @@ -156,6 +159,7 @@
156 159 character(len=PATH_MAX) :: p_e_flux_file
157 160 character(len=PATH_MAX) :: sss_file
158 161 character(len=PATH_MAX) :: sst_file
  162 + character(len=PATH_MAX) :: ice_file !DD
159 163  
160 164 REALTYPE :: wx,wy
161 165 REALTYPE :: airp
... ... @@ -166,6 +170,7 @@
166 170 REALTYPE :: const_tx,const_ty
167 171 REALTYPE :: const_swr,const_heat
168 172 REALTYPE :: const_p_e
  173 + REALTYPE :: iceatt=0.90 !DD
169 174  
170 175 REALTYPE :: es,ea,qs,qa,L,S
171 176 REALTYPE :: cdd,chd,ced
... ... @@ -262,7 +267,8 @@
262 267 momentumflux_file, &
263 268 p_e_method,const_p_e,p_e_flux_file, &
264 269 sst_method, sst_file, &
265   - sss_method, sss_file
  270 + sss_method, sss_file, &
  271 + ice_method, ice_file, iceatt !DD
266 272 !
267 273 !-----------------------------------------------------------------------
268 274 !BOC
... ... @@ -278,6 +284,16 @@
278 284 LEVEL2 'Reading meteo data from:'
279 285 LEVEL3 trim(meteo_file)
280 286  
  287 +! Sea ice !DD
  288 + select case (ice_method)
  289 + case (FROMFILE)
  290 + open(ice_unit,file=ice_file,action='read', &
  291 + status='old',err=98)
  292 + LEVEL2 'Reading ice concentration from:'
  293 + LEVEL3 trim(ice_file)
  294 + case default
  295 + end select
  296 +
281 297 else
282 298  
283 299 ! The heat fluxes
... ... @@ -336,6 +352,7 @@
336 352 cloud=0.
337 353 sss=0.
338 354 airt=0.
  355 + ice=0. !DD
339 356  
340 357 alon = deg2rad*lon
341 358 alat = deg2rad*lat
... ... @@ -358,6 +375,8 @@
358 375 stop 'init_airsea'
359 376 97 FATAL 'I could not open ',trim(sss_file)
360 377 stop 'init_airsea'
  378 +98 FATAL 'I could not open ',trim(ice_file)
  379 + stop 'init_airsea'
361 380  
362 381 end subroutine init_air_sea
363 382 !EOC
... ... @@ -395,8 +414,18 @@
395 414 !-----------------------------------------------------------------------
396 415 !BOC
397 416 if (calc_fluxes) then
  417 +! Sea ice
  418 + select case (ice_method)
  419 + case (CONSTVAL)
  420 + ice=0.50
  421 + case (FROMFILE)
  422 + call read_ice(jul,secs,ice)
  423 + case default
  424 + end select
  425 +
398 426 call flux_from_meteo(jul,secs)
399   - call short_wave_radiation(jul,secs,alon,alat)
  427 + call short_wave_radiation(jul,secs,alon,alat,ice)
  428 +
400 429 else
401 430 ! The heat fluxes
402 431 select case (heat_method)
... ... @@ -469,6 +498,7 @@
469 498 !BOC
470 499 if (calc_fluxes) then
471 500 close(meteo_unit)
  501 + if (ice_method .eq. FROMFILE) close(ice_unit) !DD
472 502 else
473 503 if (heat_method .eq. FROMFILE) close(heat_unit)
474 504 if (momentum_method .eq. FROMFILE) close(momentum_unit)
... ... @@ -737,7 +767,7 @@
737 767 ! !IROUTINE: Calculate the short--wave radiation \label{sec:swr}
738 768 !
739 769 ! !INTERFACE:
740   - subroutine short_wave_radiation(jul,secs,lon,lat,swr)
  770 + subroutine short_wave_radiation(jul,secs,lon,lat,ice,swr)
741 771 !
742 772 ! !DESCRIPTION:
743 773 ! This subroutine calculates the short--wave net radiation based on
... ... @@ -767,6 +797,7 @@
767 797 ! !INPUT PARAMETERS:
768 798 integer, intent(in) :: jul,secs
769 799 REALTYPE, intent(in) :: lon,lat
  800 + REALTYPE, intent(in) :: ice
770 801 !
771 802 ! !OUTPUT PARAMETERS:
772 803 REALTYPE, optional, intent(out) :: swr
... ... @@ -871,6 +902,11 @@
871 902  
872 903 qshort = qtot*(1.0-0.62*cloud + 0.0019*sunbet)*(1.-albedo)
873 904  
  905 +!DD Calculates the attenuation due to the presence of ice
  906 +! iceatt is a constant attenuation coefficient (default = 0.90)
  907 + qshort = qshort*(1.-iceatt*ice/100.)
  908 +
  909 +
874 910 if (present(swr)) then
875 911 swr = qshort
876 912 else
... ... @@ -956,7 +992,7 @@
956 992 call exchange_coefficients()
957 993 if (first) then
958 994 call do_calc_fluxes(heatf=h1,taux=tx1,tauy=ty1)
959   - call short_wave_radiation(jul,secs,alon,alat,swr=I1)
  995 + call short_wave_radiation(jul,secs,alon,alat,ice,swr=I1)
960 996 I2 = I1
961 997 h2 = h1
962 998 tx2 = tx1
... ... @@ -968,7 +1004,7 @@
968 1004 tx1 = tx2
969 1005 ty1 = ty2
970 1006 call do_calc_fluxes(heatf=h2,taux=tx2,tauy=ty2)
971   - call short_wave_radiation(jul,secs,alon,alat,swr=I2)
  1007 + call short_wave_radiation(jul,secs,alon,alat,ice,swr=I2)
972 1008 end if
973 1009 dt = time_diff(meteo_jul2,meteo_secs2,meteo_jul1,meteo_secs1)
974 1010 alpha(1) = (I2-I1)/dt
... ... @@ -1333,6 +1369,77 @@
1333 1369 !-----------------------------------------------------------------------
1334 1370 !BOP
1335 1371 !
  1372 +! !IROUTINE: Read sea ice, interpolate in time
  1373 +!
  1374 +! !INTERFACE:
  1375 + subroutine read_ice(jul,secs,ice)
  1376 +!
  1377 +! !DESCRIPTION:
  1378 +! For {\tt calc\_fluxes=.false.}, this routine reads sea ice
  1379 +! concentration from {\tt ice\_file}
  1380 +! and interpolates in time.
  1381 +!
  1382 +! !USES:
  1383 + IMPLICIT NONE
  1384 +!
  1385 +! !INPUT PARAMETERS:
  1386 + integer, intent(in) :: jul,secs
  1387 +!
  1388 +! !OUTPUT PARAMETERS:
  1389 + REALTYPE,intent(out) :: ice
  1390 +!
  1391 +! !REVISION HISTORY:
  1392 +! Original author(s): Karsten Bolding
  1393 +!
  1394 +! See log for airsea module
  1395 +!
  1396 +!EOP
  1397 +!
  1398 +! !LOCAL VARIABLES:
  1399 + integer :: yy,mm,dd,hh,min,ss
  1400 + REALTYPE :: t,alpha
  1401 + REALTYPE, save :: dt
  1402 + integer, save :: ice_jul1,ice_secs1
  1403 + integer, save :: ice_jul2=0,ice_secs2=0
  1404 + REALTYPE, save :: obs1(1),obs2(1)=0.
  1405 + integer :: rc
  1406 +!
  1407 +!-----------------------------------------------------------------------
  1408 +!BOC
  1409 + if (init_saved_vars) then
  1410 + ice_jul2=0
  1411 + ice_secs2=0
  1412 + obs2(1)=0.
  1413 + end if
  1414 +! This part initialise and read in new values if necessary.
  1415 + if(time_diff(ice_jul2,ice_secs2,jul,secs) .lt. 0) then
  1416 + do
  1417 + ice_jul1 = ice_jul2
  1418 + ice_secs1 = ice_secs2
  1419 + obs1 = obs2
  1420 + call read_obs(ice_unit,yy,mm,dd,hh,min,ss,1,obs2,rc)
  1421 + call julian_day(yy,mm,dd,ice_jul2)
  1422 + ice_secs2 = hh*3600 + min*60 + ss
  1423 + if(time_diff(ice_jul2,ice_secs2,jul,secs) .gt. 0) EXIT
  1424 + end do
  1425 + dt = time_diff(ice_jul2,ice_secs2,ice_jul1,ice_secs1)
  1426 + end if
  1427 +
  1428 +! Do the time interpolation
  1429 + t = time_diff(jul,secs,ice_jul1,ice_secs1)
  1430 + alpha = (obs2(1)-obs1(1))/dt
  1431 + ice = obs1(1) + t*alpha
  1432 +
  1433 + return
  1434 + end subroutine read_ice
  1435 +!EOC
  1436 +
  1437 +
  1438 +
  1439 +
  1440 +!-----------------------------------------------------------------------
  1441 +!BOP
  1442 +!
1336 1443 ! !IROUTINE: Integrate short--wave and sea surface fluxes
1337 1444 !
1338 1445 ! !INTERFACE:
... ...