source: LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 @ 2261

Last change on this file since 2261 was 2254, checked in by fhourdin, 10 years ago

Modification du calcul des flux air/mer
1) Introduction d'un facteur f_qsat_oce=0.98 devant qsat dans le calcul
de l'évaporation sur océan pour tenir compte de la moindre évaporation
de l'eau salée.
2) Introduction d'une différentiation entre z0 pour le sensible, z0h,
et le latent, z0q, imposé constant z0q=f_z0qh_oce*z0h

Modification of air/sea fluxes computation
1) Introduction of a correcting factor f_qsat_oce=0.98 on qsat
to account for the weaker evaporation of salty water.
2) Introduction of z0q=f_z0qh_oce*z0h
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 28.8 KB
RevLine 
[781]1!
2MODULE ocean_slab_mod
3!
4! This module is used for both surface ocean and sea-ice when using the slab ocean,
5! "ocean=slab".
6!
[2057]7
8  USE dimphy
9  USE indice_sol_mod
[2209]10  USE surface_data
[2057]11
[781]12  IMPLICIT NONE
[996]13  PRIVATE
[2209]14  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
[781]15
[2209]16!****************************************************************************************
17! Global saved variables
18!****************************************************************************************
19
[2057]20  INTEGER, PRIVATE, SAVE                           :: cpl_pas
21  !$OMP THREADPRIVATE(cpl_pas)
22  REAL, PRIVATE, SAVE                              :: cyang
23  !$OMP THREADPRIVATE(cyang)
24  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
25  !$OMP THREADPRIVATE(slabh)
26  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
27  !$OMP THREADPRIVATE(tslab)
[2209]28  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
29  !$OMP THREADPRIVATE(fsic)
30  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
31  !$OMP THREADPRIVATE(tice)
32  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
33  !$OMP THREADPRIVATE(seaice)
[2057]34  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
35  !$OMP THREADPRIVATE(slab_bils)
36  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
37  !$OMP THREADPRIVATE(bils_cum)
[2209]38  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
39  !$OMP THREADPRIVATE(slab_bilg)
40  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
41  !$OMP THREADPRIVATE(bilg_cum)
[2057]42
[2209]43!****************************************************************************************
44! Parameters (could be read in def file: move to slab_init)
45!****************************************************************************************
46! snow and ice physical characteristics:
47    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
48    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
49    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
50    REAL, PARAMETER :: ice_den=917.
51    REAL, PARAMETER :: sea_den=1025.
52    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
53    REAL, PARAMETER :: sno_cond=0.31*sno_den
54    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
55    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
56
57! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
58    REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
59    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
60    REAL, PARAMETER :: ice_frac_min=0.001
61    REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
62    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
63    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
64    ! below ice_thin, priority is melt lateral / grow height
65    ! ice_thin is also height of new ice
66    REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
67    ! above ice_thick, priority is melt height / grow lateral
68    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
69    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
70
71! albedo  and radiation parameters
72    REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
73    REAL, PARAMETER :: alb_sno_del=0.3  !max snow albedo = min + del
74    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
75    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
76    REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
77    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
78
79!****************************************************************************************
80
[781]81CONTAINS
82!
83!****************************************************************************************
84!
[2057]85  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
86  !, seaice_rst etc
[781]87
[2057]88    use IOIPSL
89
90    INCLUDE "iniprint.h"
91    ! For ok_xxx vars (Ekman...)
92    INCLUDE "clesphys.h"
93
94    ! Input variables
95!****************************************************************************************
96    REAL, INTENT(IN)                         :: dtime
97! Variables read from restart file
98    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
99
100! Local variables
101!****************************************************************************************
102    INTEGER                :: error
103    CHARACTER (len = 80)   :: abort_message
104    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
105
106!****************************************************************************************
107! Allocate surface fraction read from restart file
108!****************************************************************************************
[2209]109    ALLOCATE(fsic(klon), stat = error)
[2057]110    IF (error /= 0) THEN
111       abort_message='Pb allocation tmp_pctsrf_slab'
112       CALL abort_gcm(modname,abort_message,1)
113    ENDIF
[2209]114    fsic(:)=0.
115    WHERE (1.-zmasq(:)>EPSFRA)
116        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
117    END WHERE
[2057]118
119!****************************************************************************************
[2209]120! Allocate saved variables
[2057]121!****************************************************************************************
[2209]122    ALLOCATE(tslab(klon,nslay), stat=error)
123       IF (error /= 0) CALL abort_gcm &
124         (modname,'pb allocation tslab', 1)
125
[2057]126    ALLOCATE(slab_bils(klon), stat = error)
127    IF (error /= 0) THEN
128       abort_message='Pb allocation slab_bils'
129       CALL abort_gcm(modname,abort_message,1)
130    ENDIF
131    slab_bils(:) = 0.0   
132    ALLOCATE(bils_cum(klon), stat = error)
133    IF (error /= 0) THEN
134       abort_message='Pb allocation slab_bils_cum'
135       CALL abort_gcm(modname,abort_message,1)
136    ENDIF
137    bils_cum(:) = 0.0   
138
[2209]139    IF (version_ocean=='sicINT') THEN
140        ALLOCATE(slab_bilg(klon), stat = error)
141        IF (error /= 0) THEN
142           abort_message='Pb allocation slab_bilg'
143           CALL abort_gcm(modname,abort_message,1)
144        ENDIF
145        slab_bilg(:) = 0.0   
146        ALLOCATE(bilg_cum(klon), stat = error)
147        IF (error /= 0) THEN
148           abort_message='Pb allocation slab_bilg_cum'
149           CALL abort_gcm(modname,abort_message,1)
150        ENDIF
151        bilg_cum(:) = 0.0   
152        ALLOCATE(tice(klon), stat = error)
153        IF (error /= 0) THEN
154           abort_message='Pb allocation slab_tice'
155           CALL abort_gcm(modname,abort_message,1)
156        ENDIF
157        ALLOCATE(seaice(klon), stat = error)
158        IF (error /= 0) THEN
159           abort_message='Pb allocation slab_seaice'
160           CALL abort_gcm(modname,abort_message,1)
161        ENDIF
162    END IF
163
164!****************************************************************************************
165! Define some parameters
166!****************************************************************************************
[2057]167! Layer thickness
168    ALLOCATE(slabh(nslay), stat = error)
169    IF (error /= 0) THEN
170       abort_message='Pb allocation slabh'
171       CALL abort_gcm(modname,abort_message,1)
172    ENDIF
173    slabh(1)=50.
174! cyang = 1/heat capacity of top layer (rho.c.H)
175    cyang=1/(slabh(1)*4.228e+06)
176
177! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
178! pour un calcul à chaque pas de temps, cpl_pas=1
179    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
180    CALL getin('cpl_pas',cpl_pas)
181    print *,'cpl_pas',cpl_pas
[2209]182
[2057]183 END SUBROUTINE ocean_slab_init
184!
185!****************************************************************************************
186!
187  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
188
[996]189    USE limit_read_mod
[1785]190
[996]191!    INCLUDE "clesphys.h"
[781]192
[996]193! Arguments
[781]194!****************************************************************************************
[996]195    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
196    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
197    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
[2057]198    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
[996]199    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
[781]200
201! Local variables
202!****************************************************************************************
203
[2057]204    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
205       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
[996]206    ELSE
[2209]207       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
208       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
[2057]209       is_modified=.TRUE.
[996]210    END IF
[781]211
[996]212  END SUBROUTINE ocean_slab_frac
[781]213!
214!****************************************************************************************
215!
216  SUBROUTINE ocean_slab_noice( &
[996]217       itime, dtime, jour, knon, knindex, &
[2254]218       p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, &
[1067]219       AcoefH, AcoefQ, BcoefH, BcoefQ, &
220       AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]221       ps, u1, v1, gustiness, tsurf_in, &
[2209]222       radsol, snow, &
[1067]223       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[2057]224       tsurf_new, dflux_s, dflux_l, qflux)
[1067]225   
226    USE calcul_fluxs_mod
[1785]227
[781]228    INCLUDE "iniprint.h"
[2254]229    INCLUDE "clesphys.h"
[781]230
231! Input arguments
232!****************************************************************************************
[996]233    INTEGER, INTENT(IN)                  :: itime
234    INTEGER, INTENT(IN)                  :: jour
[781]235    INTEGER, INTENT(IN)                  :: knon
236    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
237    REAL, INTENT(IN)                     :: dtime
238    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
[2254]239    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
[781]240    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
241    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
[1067]242    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
243    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
[781]244    REAL, DIMENSION(klon), INTENT(IN)    :: ps
[2240]245    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
[996]246    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
[2209]247    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
[781]248
249! In/Output arguments
250!****************************************************************************************
251    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
252   
253! Output arguments
254!****************************************************************************************
255    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
256    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
[1067]257    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
[781]258    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
259    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
[2057]260    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
[781]261
262! Local variables
263!****************************************************************************************
[2057]264    INTEGER               :: i,ki
[2209]265    ! surface fluxes
[996]266    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
[2209]267    ! flux correction
268    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
269    ! surface wind stress
[1067]270    REAL, DIMENSION(klon) :: u0, v0
271    REAL, DIMENSION(klon) :: u1_lay, v1_lay
[2209]272    ! ice creation
273    REAL                  :: e_freeze, h_new, dfsic
[781]274
275!****************************************************************************************
[996]276! 1) Flux calculation
277!
278!****************************************************************************************
279    cal(:)      = 0.
280    beta(:)     = 1.
281    dif_grnd(:) = 0.
[781]282   
[1067]283! Suppose zero surface speed
284    u0(:)=0.0
285    v0(:)=0.0
286    u1_lay(:) = u1(:) - u0(:)
287    v1_lay(:) = v1(:) - v0(:)
288
[781]289    CALL calcul_fluxs(knon, is_oce, dtime, &
[2254]290         tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
[781]291         precip_rain, precip_snow, snow, qsurf,  &
[2240]292         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
[2254]293         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]294         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
295
[1067]296! - Flux calculation at first modele level for U and V
297    CALL calcul_flux_wind(knon, dtime, &
[2240]298         u0, v0, u1, v1, gustiness, cdragm, &
[1067]299         AcoefU, AcoefV, BcoefU, BcoefV, &
300         p1lay, temp_air, &
301         flux_u1, flux_v1) 
302
[2057]303! Accumulate total fluxes locally
304    slab_bils(:)=0.
305    DO i=1,knon
306        ki=knindex(i)
[2209]307        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
308                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
[2057]309        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
310! Also taux, tauy, saved vars...
311    END DO
312
[781]313!****************************************************************************************
[2057]314! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
[781]315!
316!****************************************************************************************
[2057]317    lmt_bils(:)=0.
[2209]318    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
319    ! lmt_bils and diff_sst,siv saved by limit_slab
320    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
[2057]321    ! qflux = total QFlux correction (in W/m2)
[781]322
323!****************************************************************************************
[996]324! 3) Recalculate new temperature
[781]325!
[2057]326!***********************************************o*****************************************
327    tsurf_new=tsurf_in
328    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
329        ! Compute transport
330        ! Add read QFlux and SST tendency
331        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
332        ! Add cumulated surface fluxes
333        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
334        ! Update surface temperature
335        SELECT CASE(version_ocean)
336        CASE('sicNO')
337            DO i=1,knon
338                ki=knindex(i)
339                tsurf_new(i)=tslab(ki,1)
340            END DO
[2209]341        CASE('sicOBS') ! check for sea ice or tslab below freezing
[2057]342            DO i=1,knon
343                ki=knindex(i)
[2209]344                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
[2057]345                    tslab(ki,1)=t_freeze
346                END IF
[2209]347                tsurf_new(i)=tslab(ki,1)
[2057]348            END DO
349        CASE('sicINT')
350            DO i=1,knon
351                ki=knindex(i)
[2209]352                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
353                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
354                        ! quantity of new ice formed
355                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
356                        ! new ice
357                        tice(ki)=t_freeze
358                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
359                        IF (fsic(ki).GT.ice_frac_min) THEN
360                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
361                            tslab(ki,1)=t_freeze
362                        ELSE
363                            fsic(ki)=0.
364                        END IF
365                        tsurf_new(i)=t_freeze
366                    ELSE
[2057]367                        tsurf_new(i)=tslab(ki,1)
[2209]368                    END IF
369                ELSE ! ice present
[2057]370                    tsurf_new(i)=t_freeze
[2209]371                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
372                        ! quantity of new ice formed over open ocean
373                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
374                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
375                        ! new ice height and fraction
376                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
377                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
378                        h_new=MIN(e_freeze/dfsic,h_ice_max)
379                        ! update tslab to freezing over open ocean only
380                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
381                        ! update sea ice
382                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
383                                   /(dfsic+fsic(ki))
384                        fsic(ki)=fsic(ki)+dfsic
385                        ! update snow?
386                    END IF !freezing
387                END IF ! sea ice present
[2057]388            END DO
389        END SELECT
390        bils_cum(:)=0.0! clear cumulated fluxes
391    END IF ! coupling time
392  END SUBROUTINE ocean_slab_noice
393!
[781]394!****************************************************************************************
[2209]395
396  SUBROUTINE ocean_slab_ice(   &
397       itime, dtime, jour, knon, knindex, &
398       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
399       AcoefH, AcoefQ, BcoefH, BcoefQ, &
400       AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]401       ps, u1, v1, gustiness, &
[2209]402       radsol, snow, qsurf, qsol, agesno, &
403       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
404       tsurf_new, dflux_s, dflux_l, swnet)
405
406   USE calcul_fluxs_mod
407
408   INCLUDE "YOMCST.h"
[2254]409   INCLUDE "clesphys.h"
[2209]410
411! Input arguments
[2057]412!****************************************************************************************
[2209]413    INTEGER, INTENT(IN)                  :: itime, jour, knon
414    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
415    REAL, INTENT(IN)                     :: dtime
416    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
417    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
418    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
419    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
420    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
421    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
422    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
423    REAL, DIMENSION(klon), INTENT(IN)    :: ps
[2240]424    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
[2209]425    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
426
427! In/Output arguments
428!****************************************************************************************
429    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
430    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
431    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
432
433! Output arguments
434!****************************************************************************************
435    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
436    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
437    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
438    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
439    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
440    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
441    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l
442
443! Local variables
444!****************************************************************************************
445    INTEGER               :: i,ki
446    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
447    REAL, DIMENSION(klon) :: u0, v0
448    REAL, DIMENSION(klon) :: u1_lay, v1_lay
449    ! intermediate heat fluxes:
450    REAL                  :: f_cond, f_swpen
451    ! for snow/ice albedo:
452    REAL                  :: alb_snow, alb_ice, alb_pond
453    REAL                  :: frac_snow, frac_ice, frac_pond
454    ! for ice melt / freeze
455    REAL                  :: e_melt, snow_evap, h_test
456    ! dhsic, dfsic change in ice mass, fraction.
457    REAL                  :: dhsic, dfsic, frac_mf
458
459!****************************************************************************************
[2057]460! 1) Flux calculation
461!****************************************************************************************
[2209]462! Suppose zero surface speed
463    u0(:)=0.0
464    v0(:)=0.0
465    u1_lay(:) = u1(:) - u0(:)
466    v1_lay(:) = v1(:) - v0(:)
467
468! set beta, cal, compute conduction fluxes inside ice/snow
469    slab_bilg(:)=0.
470    dif_grnd(:)=0.
471    beta(:) = 1.
472    DO i=1,knon
473    ki=knindex(i)
474        IF (snow(i).GT.snow_min) THEN
475            ! snow-layer heat capacity
476            cal(i)=2.*RCPD/(snow(i)*ice_cap)
477            ! snow conductive flux
478            f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
479            ! all shortwave flux absorbed
480            f_swpen=0.
481            ! bottom flux (ice conduction)
482            slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
483            ! update ice temperature
484            tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
485                     *(slab_bilg(ki)+f_cond)*dtime
486       ELSE ! bare ice
487            ! ice-layer heat capacity
488            cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
489            ! conductive flux
490            f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
491            ! penetrative shortwave flux...
492            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
493            slab_bilg(ki)=f_swpen-f_cond
494        END IF
495        radsol(i)=radsol(i)+f_cond-f_swpen
496    END DO
497    ! weight fluxes to ocean by sea ice fraction
498    slab_bilg(:)=slab_bilg(:)*fsic(:)
499
[2057]500! calcul_fluxs (sens, lat etc)
[2209]501    CALL calcul_fluxs(knon, is_sic, dtime, &
[2254]502        tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
[2209]503        precip_rain, precip_snow, snow, qsurf,  &
[2240]504        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
[2254]505        f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[2209]506        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
507    DO i=1,knon
508        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
509    END DO
510
[2057]511! calcul_flux_wind
[2209]512    CALL calcul_flux_wind(knon, dtime, &
[2240]513         u0, v0, u1, v1, gustiness, cdragm, &
[2209]514         AcoefU, AcoefV, BcoefU, BcoefV, &
515         p1lay, temp_air, &
516         flux_u1, flux_v1)
[781]517
[2057]518!****************************************************************************************
[2209]519! 2) Update snow and ice surface
[2057]520!****************************************************************************************
[2209]521! snow precip
522    DO i=1,knon
523        ki=knindex(i)
524        IF (precip_snow(i) > 0.) THEN
525            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
526        END IF
527! snow and ice sublimation
528        IF (evap(i) > 0.) THEN
529           snow_evap = MIN (snow(i) / dtime, evap(i))
530           snow(i) = snow(i) - snow_evap * dtime
531           snow(i) = MAX(0.0, snow(i))
532           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
533        ENDIF
534! Melt / Freeze from above if Tsurf>0
535        IF (tsurf_new(i).GT.t_melt) THEN
536            ! energy available for melting snow (in kg/m2 of snow)
537            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
538               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
539            ! remove snow
540            IF (snow(i).GT.e_melt) THEN
541                snow(i)=snow(i)-e_melt
542                tsurf_new(i)=t_melt
543            ELSE ! all snow is melted
544                ! add remaining heat flux to ice
545                e_melt=e_melt-snow(i)
546                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
547                tsurf_new(i)=tice(ki)
548            END IF
549        END IF
550! melt ice from above if Tice>0
551        IF (tice(ki).GT.t_melt) THEN
552            ! quantity of ice melted (kg/m2)
553            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
554             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
555            ! melt from above, height only
556            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
557            e_melt=e_melt-dhsic
558            IF (e_melt.GT.0) THEN
559            ! lateral melt if ice too thin
560            dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
561            ! if all melted add remaining heat to ocean
562            e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
563            slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
564            ! update height and fraction
565            fsic(ki)=fsic(ki)-dfsic
566            END IF
567            seaice(ki)=seaice(ki)-dhsic
568            ! surface temperature at melting point
569            tice(ki)=t_melt
570            tsurf_new(i)=t_melt
571        END IF
572        ! convert snow to ice if below floating line
573        h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den
574        IF (h_test.GT.0.) THEN !snow under water
575            ! extra snow converted to ice (with added frozen sea water)
576            dhsic=h_test/(sea_den-ice_den+sno_den)
577            seaice(ki)=seaice(ki)+dhsic
578            snow(i)=snow(i)-dhsic*sno_den/ice_den
579            ! available energy (freeze sea water + bring to tice)
580            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
581                   ice_cap/2.*(t_freeze-tice(ki)))
582            ! update ice temperature
583            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
584        END IF
585    END DO
586
[2057]587! New albedo
[2209]588    DO i=1,knon
589        ki=knindex(i)
590       ! snow albedo: update snow age
591        IF (snow(i).GT.0.0001) THEN
592             agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
593                         * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.)
594        ELSE
595            agesno(i)=0.0
596        END IF
597        ! snow albedo
598        alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
599        ! ice albedo (varies with ice tkickness and temp)
600        alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
601        IF (tice(ki).GT.t_freeze-0.01) THEN
602            alb_ice=MIN(alb_ice,alb_ice_wet)
603        ELSE
604            alb_ice=MIN(alb_ice,alb_ice_dry)
605        END IF
606        ! pond albedo
607        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
608        ! pond fraction
609        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
610        ! snow fraction
611        frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
612        ! ice fraction
613        frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
614        ! total albedo
615        alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
616    END DO
617    alb2_new(:) = alb1_new(:)
[2057]618
619!****************************************************************************************
[2209]620! 3) Recalculate new ocean temperature (add fluxes below ice)
[2057]621!    Melt / freeze from below
622!***********************************************o*****************************************
[2209]623    !cumul fluxes
624    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
625    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
626        ! Add cumulated surface fluxes
627        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
628        DO i=1,knon
629            ki=knindex(i)
630            ! split lateral/top melt-freeze
631            frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
632            IF (tslab(ki,1).LE.t_freeze) THEN
633               ! ****** Form new ice from below *******
634               ! quantity of new ice
635                e_melt=(t_freeze-tslab(ki,1))/cyang &
636                       /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
637               ! first increase height to h_thin
638               dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
639               seaice(ki)=dhsic+seaice(ki)
640               e_melt=e_melt-fsic(ki)*dhsic
641               IF (e_melt.GT.0.) THEN
642               ! frac_mf fraction used for lateral increase
643               dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
644               fsic(ki)=fsic(ki)+dfsic
645               e_melt=e_melt-dfsic*seaice(ki)
646               ! rest used to increase height
647               seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
648               END IF
649               tslab(ki,1)=t_freeze
650           ELSE ! slab temperature above freezing
651               ! ****** melt ice from below *******
652               ! quantity of melted ice
653               e_melt=(tslab(ki,1)-t_freeze)/cyang &
654                       /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
655               ! first decrease height to h_thick
656               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
657               seaice(ki)=seaice(ki)-dhsic
658               e_melt=e_melt-fsic(ki)*dhsic
659               IF (e_melt.GT.0) THEN
660               ! frac_mf fraction used for height decrease
661               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
662               seaice(ki)=seaice(ki)-dhsic
663               e_melt=e_melt-fsic(ki)*dhsic
664               ! rest used to decrease fraction (up to 0!)
665               dfsic=MIN(fsic(ki),e_melt/seaice(ki))
666               ! keep remaining in ocean
667               e_melt=e_melt-dfsic*seaice(ki)
668               END IF
669               tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang
670               fsic(ki)=fsic(ki)-dfsic
671           END IF
672        END DO
673        bilg_cum(:)=0.
674    END IF ! coupling time
675   
676    !tests fraction glace
677    WHERE (fsic.LT.ice_frac_min)
678        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
679        tice=t_melt
680        fsic=0.
681        seaice=0.
682    END WHERE
[2057]683
[2209]684  END SUBROUTINE ocean_slab_ice
[781]685!
686!****************************************************************************************
687!
[2057]688  SUBROUTINE ocean_slab_final
689
690!****************************************************************************************
691! Deallocate module variables
692!****************************************************************************************
693    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
[2209]694    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
695    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
696    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
697    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
698    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
699    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
[2057]700
701  END SUBROUTINE ocean_slab_final
702
[781]703END MODULE ocean_slab_mod
Note: See TracBrowser for help on using the repository browser.