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

Last change on this file since 2989 was 2656, checked in by Ehouarn Millour, 8 years ago

Making the slab work:

  • added a slab_heat_transp_mod module for horizontal diffusion and Ekman transport
  • added storage and output of relevent variables in phyredem, phyetat0, phy_output_ctrlout_mod, phys_output_write_mod
  • moved nslay (number of slab layers) out of dimphy into ocean_slab_mod.

FC

  • 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: 37.5 KB
Line 
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!
7
8  USE dimphy
9  USE indice_sol_mod
10  USE surface_data
11  USE mod_grid_phy_lmdz, ONLY: klon_glo
12  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
13
14  IMPLICIT NONE
15  PRIVATE
16  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
17
18!***********************************************************************************
19! Global saved variables
20!***********************************************************************************
21  ! number of slab vertical layers
22  INTEGER, PUBLIC, SAVE :: nslay
23  !$OMP THREADPRIVATE(nslay)
24  ! timestep for coupling (update slab temperature) in timesteps
25  INTEGER, PRIVATE, SAVE                           :: cpl_pas
26  !$OMP THREADPRIVATE(cpl_pas)
27  ! cyang = 1/heat capacity of top layer (rho.c.H)
28  REAL, PRIVATE, SAVE                              :: cyang
29  !$OMP THREADPRIVATE(cyang)
30  ! depth of slab layers (1 or 2)
31  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
32  !$OMP THREADPRIVATE(slabh)
33  ! slab temperature
34  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
35  !$OMP THREADPRIVATE(tslab)
36  ! heat flux convergence due to Ekman
37  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
38  !$OMP THREADPRIVATE(dt_ekman)
39  ! heat flux convergence due to horiz diffusion
40  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
41  !$OMP THREADPRIVATE(dt_hdiff)
42  ! fraction of ocean covered by sea ice (sic / (oce+sic))
43  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
44  !$OMP THREADPRIVATE(fsic)
45  ! temperature of the sea ice
46  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
47  !$OMP THREADPRIVATE(tice)
48  ! sea ice thickness, in kg/m2
49  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
50  !$OMP THREADPRIVATE(seaice)
51  ! net surface heat flux, weighted by open ocean fraction
52  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
53  !$OMP THREADPRIVATE(slab_bils)
54  ! slab_bils accumulated over cpl_pas timesteps
55  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
56  !$OMP THREADPRIVATE(bils_cum)
57  ! net heat flux into the ocean below the ice : conduction + solar radiation
58  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
59  !$OMP THREADPRIVATE(slab_bilg)
60  ! slab_bilg over cpl_pas timesteps
61  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
62  !$OMP THREADPRIVATE(bilg_cum)
63  ! wind stress saved over cpl_pas timesteps
64  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
65  !$OMP THREADPRIVATE(taux_cum)
66  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
67  !$OMP THREADPRIVATE(tauy_cum)
68
69!***********************************************************************************
70! Parameters (could be read in def file: move to slab_init)
71!***********************************************************************************
72! snow and ice physical characteristics:
73    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
74    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
75    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
76    REAL, PARAMETER :: ice_den=917. ! ice density
77    REAL, PARAMETER :: sea_den=1025. ! sea water density
78    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
79    REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
80    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
81    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, snow and ice
82    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
83
84! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
85    REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
86    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
87    REAL, PARAMETER :: ice_frac_min=0.001
88    REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
89    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
90    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
91    ! below ice_thin, priority is melt lateral / grow height
92    ! ice_thin is also height of new ice
93    REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
94    ! above ice_thick, priority is melt height / grow lateral
95    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
96    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
97
98! albedo  and radiation parameters
99    REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
100    REAL, PARAMETER :: alb_sno_del=0.3  !max snow albedo = min + del
101    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
102    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
103    REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
104    ! ice (no snow)
105    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
106
107! horizontal transport
108   LOGICAL, PUBLIC, SAVE :: slab_hdiff
109   !$OMP THREADPRIVATE(slab_hdiff)
110   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
111   !$OMP THREADPRIVATE(coef_hdiff)
112   INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
113   !$OMP THREADPRIVATE(slab_ekman)
114   !$OMP THREADPRIVATE(slab_cadj)
115
116!***********************************************************************************
117
118CONTAINS
119!
120!***********************************************************************************
121!
122  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
123  !, seaice_rst etc
124
125    USE ioipsl_getin_p_mod, ONLY : getin_p
126    USE mod_phys_lmdz_transfert_para, ONLY : gather
127    USE slab_heat_transp_mod, ONLY : ini_slab_transp
128
129    ! Input variables
130!***********************************************************************************
131    REAL, INTENT(IN)                         :: dtime
132! Variables read from restart file
133    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
134    ! surface fractions from start file
135
136! Local variables
137!************************************************************************************
138    INTEGER                :: error
139    REAL, DIMENSION(klon_glo) :: zmasq_glo
140    CHARACTER (len = 80)   :: abort_message
141    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
142
143!***********************************************************************************
144! Define some parameters
145!***********************************************************************************
146! Number of slab layers
147    nslay=2
148    CALL getin_p('slab_layers',nslay)
149    print *,'number of slab layers : ',nslay
150! Layer thickness
151    ALLOCATE(slabh(nslay), stat = error)
152    IF (error /= 0) THEN
153       abort_message='Pb allocation slabh'
154       CALL abort_physic(modname,abort_message,1)
155    ENDIF
156    slabh(1)=50.
157    IF (nslay.GT.1) THEN
158        slabh(2)=150.
159    END IF
160
161! cyang = 1/heat capacity of top layer (rho.c.H)
162    cyang=1/(slabh(1)*sea_den*sea_cap)
163
164! cpl_pas  coupling period (update of tslab and ice fraction)
165! pour un calcul a chaque pas de temps, cpl_pas=1
166    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
167    CALL getin_p('cpl_pas',cpl_pas)
168    print *,'cpl_pas',cpl_pas
169
170! Horizontal diffusion
171    slab_hdiff=.FALSE.
172    CALL getin_p('slab_hdiff',slab_hdiff)
173    coef_hdiff=25000.
174    CALL getin_p('coef_hdiff',coef_hdiff)
175! Ekman transport
176    slab_ekman=0
177    CALL getin_p('slab_ekman',slab_ekman)
178! Convective adjustment
179    IF (nslay.EQ.1) THEN
180        slab_cadj=0
181    ELSE
182        slab_cadj=1
183    END IF
184    CALL getin_p('slab_cadj',slab_cadj)
185
186!************************************************************************************
187! Allocate surface fraction read from restart file
188!************************************************************************************
189    ALLOCATE(fsic(klon), stat = error)
190    IF (error /= 0) THEN
191       abort_message='Pb allocation tmp_pctsrf_slab'
192       CALL abort_physic(modname,abort_message,1)
193    ENDIF
194    fsic(:)=0.
195    !zmasq = continent fraction
196    WHERE (1.-zmasq(:)>EPSFRA)
197        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
198    END WHERE
199
200!************************************************************************************
201! Allocate saved fields
202!************************************************************************************
203    ALLOCATE(tslab(klon,nslay), stat=error)
204       IF (error /= 0) CALL abort_physic &
205         (modname,'pb allocation tslab', 1)
206
207    ALLOCATE(slab_bils(klon), stat = error)
208    IF (error /= 0) THEN
209       abort_message='Pb allocation slab_bils'
210       CALL abort_physic(modname,abort_message,1)
211    ENDIF
212    slab_bils(:) = 0.0   
213    ALLOCATE(bils_cum(klon), stat = error)
214    IF (error /= 0) THEN
215       abort_message='Pb allocation slab_bils_cum'
216       CALL abort_physic(modname,abort_message,1)
217    ENDIF
218    bils_cum(:) = 0.0   
219
220    IF (version_ocean=='sicINT') THEN ! interactive sea ice
221        ALLOCATE(slab_bilg(klon), stat = error)
222        IF (error /= 0) THEN
223           abort_message='Pb allocation slab_bilg'
224           CALL abort_physic(modname,abort_message,1)
225        ENDIF
226        slab_bilg(:) = 0.0   
227        ALLOCATE(bilg_cum(klon), stat = error)
228        IF (error /= 0) THEN
229           abort_message='Pb allocation slab_bilg_cum'
230           CALL abort_physic(modname,abort_message,1)
231        ENDIF
232        bilg_cum(:) = 0.0   
233        ALLOCATE(tice(klon), stat = error)
234        IF (error /= 0) THEN
235           abort_message='Pb allocation slab_tice'
236           CALL abort_physic(modname,abort_message,1)
237        ENDIF
238        ALLOCATE(seaice(klon), stat = error)
239        IF (error /= 0) THEN
240           abort_message='Pb allocation slab_seaice'
241           CALL abort_physic(modname,abort_message,1)
242        ENDIF
243    END IF
244
245    IF (slab_hdiff) THEN !horizontal diffusion
246        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
247        IF (error /= 0) THEN
248           abort_message='Pb allocation dt_hdiff'
249           CALL abort_physic(modname,abort_message,1)
250        ENDIF
251        dt_hdiff(:,:) = 0.0   
252    ENDIF
253
254    IF (slab_ekman.GT.0) THEN ! ekman transport
255        ALLOCATE(dt_ekman(klon,nslay), stat = error)
256        IF (error /= 0) THEN
257           abort_message='Pb allocation dt_ekman'
258           CALL abort_physic(modname,abort_message,1)
259        ENDIF
260        dt_ekman(:,:) = 0.0   
261        ALLOCATE(taux_cum(klon), stat = error)
262        IF (error /= 0) THEN
263           abort_message='Pb allocation taux_cum'
264           CALL abort_physic(modname,abort_message,1)
265        ENDIF
266        taux_cum(:) = 0.0   
267        ALLOCATE(tauy_cum(klon), stat = error)
268        IF (error /= 0) THEN
269           abort_message='Pb allocation tauy_cum'
270           CALL abort_physic(modname,abort_message,1)
271        ENDIF
272        tauy_cum(:) = 0.0   
273    ENDIF
274
275! Initialize transport
276    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
277      CALL gather(zmasq,zmasq_glo)
278!$OMP MASTER  ! Only master thread
279      IF (is_mpi_root) THEN
280          CALL ini_slab_transp(zmasq_glo)
281      END IF
282!$OMP END MASTER
283    END IF
284
285 END SUBROUTINE ocean_slab_init
286!
287!***********************************************************************************
288!
289  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
290
291! this routine sends back the sea ice and ocean fraction to the main physics
292! routine. Called only with interactive sea ice
293
294! Arguments
295!************************************************************************************
296    INTEGER, INTENT(IN)                        :: itime   ! current timestep
297    INTEGER, INTENT(IN)                        :: jour    !  day in year (not
298    REAL   , INTENT(IN)                        :: dtime   ! physics timestep (s)
299    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
300    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is
301                                                         ! modified at this time step
302
303       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
304       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
305       is_modified=.TRUE.
306
307  END SUBROUTINE ocean_slab_frac
308!
309!************************************************************************************
310!
311  SUBROUTINE ocean_slab_noice( &
312       itime, dtime, jour, knon, knindex, &
313       p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, &
314       AcoefH, AcoefQ, BcoefH, BcoefQ, &
315       AcoefU, AcoefV, BcoefU, BcoefV, &
316       ps, u1, v1, gustiness, tsurf_in, &
317       radsol, snow, &
318       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
319       tsurf_new, dflux_s, dflux_l, qflux)
320   
321    USE calcul_fluxs_mod
322    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
323    USE mod_phys_lmdz_para
324
325    INCLUDE "clesphys.h"
326
327! This routine
328! (1) computes surface turbulent fluxes over points with some open ocean   
329! (2) reads additional Q-flux (everywhere)
330! (3) computes horizontal transport (diffusion & Ekman)
331! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
332
333! Note :
334! klon total number of points
335! knon number of points with open ocean (varies with time)
336! knindex gives position of the knon points within klon.
337! In general, local saved variables have klon values
338! variables exchanged with PBL module have knon.
339
340! Input arguments
341!***********************************************************************************
342    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
343    INTEGER, INTENT(IN)                  :: jour  ! day in year (for Q-Flux)
344    INTEGER, INTENT(IN)                  :: knon  ! number of points
345    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
346    REAL, INTENT(IN) :: dtime  ! timestep (s)
347    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
348    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
349    ! drag coefficients
350    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
351    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum ! near surface T, q
352    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
353    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
354    ! exchange coefficients for boundary layer scheme
355    REAL, DIMENSION(klon), INTENT(IN)    :: ps  ! surface pressure
356    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness ! surface wind
357    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
358    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
359
360! In/Output arguments
361!************************************************************************************
362    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
363   
364! Output arguments
365!************************************************************************************
366    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
367    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
368    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
369    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
370    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
371    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
372
373! Local variables
374!************************************************************************************
375    INTEGER               :: i,ki,k
376    REAL                  :: t_cadj
377    !  for surface heat fluxes
378    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
379    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes
380    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
381    ! for surface wind stress
382    REAL, DIMENSION(klon) :: u0, v0
383    REAL, DIMENSION(klon) :: u1_lay, v1_lay
384    ! for new ice creation
385    REAL                  :: e_freeze, h_new, dfsic
386    ! horizontal diffusion and Ekman local vars
387    ! dimension = global domain (klon_glo) instead of // subdomains
388    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
389    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
390    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
391    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
392    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
393
394!****************************************************************************************
395! 1) Surface fluxes calculation
396!   
397!****************************************************************************************
398    cal(:)      = 0. ! infinite thermal inertia
399    beta(:)     = 1. ! wet surface
400    dif_grnd(:) = 0. ! no diffusion into ground
401   
402! Suppose zero surface speed
403    u0(:)=0.0
404    v0(:)=0.0
405    u1_lay(:) = u1(:) - u0(:)
406    v1_lay(:) = v1(:) - v0(:)
407
408! Compute latent & sensible fluxes
409    CALL calcul_fluxs(knon, is_oce, dtime, &
410         tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
411         precip_rain, precip_snow, snow, qsurf,  &
412         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
413         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
414         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
415
416! save total cumulated heat fluxes locally
417! radiative + turbulent + melt of falling snow
418    slab_bils(:)=0.
419    DO i=1,knon
420        ki=knindex(i)
421        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
422                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
423        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
424    END DO
425
426!  Compute surface wind stress
427    CALL calcul_flux_wind(knon, dtime, &
428         u0, v0, u1, v1, gustiness, cdragm, &
429         AcoefU, AcoefV, BcoefU, BcoefV, &
430         p1lay, temp_air, &
431         flux_u1, flux_v1) 
432
433! save cumulated wind stress
434    IF (slab_ekman.GT.0) THEN
435      DO i=1,knon
436          ki=knindex(i)
437          taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
438          tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
439      END DO
440    ENDIF
441
442!****************************************************************************************
443! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
444!
445!****************************************************************************************
446    lmt_bils(:)=0.
447    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
448    ! lmt_bils and diff_sst,siv saved by limit_slab
449    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
450    ! qflux = total QFlux correction (in W/m2)
451
452!****************************************************************************************
453! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
454!    Bring to freezing temp and make sea ice if necessary
455
456!***********************************************o*****************************************
457    tsurf_new=tsurf_in
458    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
459! ***********************************
460!  Horizontal transport
461! ***********************************
462      IF (slab_ekman.GT.0) THEN
463          ! copy wind stress to global var
464          CALL gather(taux_cum,taux_glo)
465          CALL gather(tauy_cum,tauy_glo)
466      END IF
467
468      IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
469        CALL gather(tslab,tslab_glo)
470      ! Compute horiz transport on one process only
471!$OMP MASTER  ! Only master thread
472        IF (is_mpi_root) THEN ! Only master processus         
473          IF (slab_hdiff) THEN
474              dt_hdiff_glo(:,:)=0.
475          END IF
476          IF (slab_ekman.GT.0) THEN
477              dt_ekman_glo(:,:)=0.
478          END IF
479          DO i=1,cpl_pas ! time splitting for numerical stability
480            IF (slab_ekman.GT.0) THEN
481              SELECT CASE (slab_ekman)
482                CASE (1)
483                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
484                CASE (2)
485                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
486                CASE DEFAULT
487                  dt_ekman_tmp(:,:)=0.
488              END SELECT
489              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
490              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1   
491              DO k=1,nslay
492                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
493              ENDDO
494              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
495            ENDIF
496            IF (slab_hdiff) THEN ! horizontal diffusion
497              ! laplacian of slab T
498              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
499              ! multiply by diff coef and normalize to 50m slab equivalent
500              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
501              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
502              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
503            END IF
504          END DO ! time splitting
505          IF (slab_hdiff) THEN
506            !dt_hdiff_glo saved in W/m2
507            DO k=1,nslay
508              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
509            END DO
510          END IF
511          IF (slab_ekman.GT.0) THEN
512            ! dt_ekman_glo saved in W/m2
513            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
514          END IF
515        END IF ! is_mpi_root
516!$OMP END MASTER
517!$OMP BARRIER
518        ! Send new fields back to all processes
519        CALL Scatter(tslab_glo,tslab)
520        IF (slab_hdiff) THEN
521          CALL Scatter(dt_hdiff_glo,dt_hdiff)
522        END IF
523        IF (slab_ekman.GT.0) THEN
524          CALL Scatter(dt_ekman_glo,dt_ekman)
525          ! clear wind stress
526          taux_cum(:)=0.
527          tauy_cum(:)=0.
528        END IF
529      ENDIF ! transport
530
531! ***********************************
532! Other heat fluxes
533! ***********************************
534      ! Add read QFlux
535      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
536      ! Add cumulated surface fluxes
537      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
538      ! Convective adjustment if 2 layers
539      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
540        DO i=1,klon
541          IF (tslab(i,2).GT.tslab(i,1)) THEN
542            ! mean (mass-weighted) temperature
543            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
544            tslab(i,1)=t_cadj
545            tslab(i,2)=t_cadj
546          END IF
547        END DO
548      END IF
549! ***********************************
550! Update surface temperature and ice
551! ***********************************
552      SELECT CASE(version_ocean)
553      CASE('sicNO') ! no sea ice even below freezing !
554          DO i=1,knon
555              ki=knindex(i)
556              tsurf_new(i)=tslab(ki,1)
557          END DO
558      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
559        ! tslab cannot be below freezing, or above it if there is sea ice
560          DO i=1,knon
561              ki=knindex(i)
562              IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
563                  tslab(ki,1)=t_freeze
564              END IF
565              tsurf_new(i)=tslab(ki,1)
566          END DO
567      CASE('sicINT') ! interactive sea ice
568          DO i=1,knon
569              ki=knindex(i)
570              IF (fsic(ki).LT.epsfra) THEN ! Free of ice
571                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
572                      ! quantity of new ice formed
573                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
574                      ! new ice
575                      tice(ki)=t_freeze
576                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
577                      IF (fsic(ki).GT.ice_frac_min) THEN
578                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
579                          tslab(ki,1)=t_freeze
580                      ELSE
581                          fsic(ki)=0.
582                      END IF
583                      tsurf_new(i)=t_freeze
584                  ELSE
585                      tsurf_new(i)=tslab(ki,1)
586                  END IF
587              ELSE ! ice present
588                  tsurf_new(i)=t_freeze
589                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
590                      ! quantity of new ice formed over open ocean
591                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
592                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
593                      ! new ice height and fraction
594                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
595                      dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
596                      h_new=MIN(e_freeze/dfsic,h_ice_max)
597                      ! update tslab to freezing over open ocean only
598                      tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
599                      ! update sea ice
600                      seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
601                                 /(dfsic+fsic(ki))
602                      fsic(ki)=fsic(ki)+dfsic
603                      ! update snow?
604                  END IF ! tslab below freezing
605              END IF ! sea ice present
606          END DO
607      END SELECT
608      bils_cum(:)=0.0! clear cumulated fluxes
609    END IF ! coupling time
610  END SUBROUTINE ocean_slab_noice
611!
612!****************************************************************************************
613
614  SUBROUTINE ocean_slab_ice(   &
615       itime, dtime, jour, knon, knindex, &
616       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
617       AcoefH, AcoefQ, BcoefH, BcoefQ, &
618       AcoefU, AcoefV, BcoefU, BcoefV, &
619       ps, u1, v1, gustiness, &
620       radsol, snow, qsurf, qsol, agesno, &
621       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
622       tsurf_new, dflux_s, dflux_l, swnet)
623
624   USE calcul_fluxs_mod
625
626   INCLUDE "YOMCST.h"
627   INCLUDE "clesphys.h"
628
629! Input arguments
630!****************************************************************************************
631    INTEGER, INTENT(IN)                  :: itime, jour, knon
632    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
633    REAL, INTENT(IN)                     :: dtime
634    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
635    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
636    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
637    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
638    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
639    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
640    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
641    REAL, DIMENSION(klon), INTENT(IN)    :: ps
642    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
643    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
644
645! In/Output arguments
646!****************************************************************************************
647    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
648    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
649    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
650
651! Output arguments
652!****************************************************************************************
653    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
654    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
655    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
656    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
657    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
658    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
659    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l
660
661! Local variables
662!****************************************************************************************
663    INTEGER               :: i,ki
664    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
665    REAL, DIMENSION(klon) :: u0, v0
666    REAL, DIMENSION(klon) :: u1_lay, v1_lay
667    ! intermediate heat fluxes:
668    REAL                  :: f_cond, f_swpen
669    ! for snow/ice albedo:
670    REAL                  :: alb_snow, alb_ice, alb_pond
671    REAL                  :: frac_snow, frac_ice, frac_pond
672    ! for ice melt / freeze
673    REAL                  :: e_melt, snow_evap, h_test
674    ! dhsic, dfsic change in ice mass, fraction.
675    REAL                  :: dhsic, dfsic, frac_mf
676
677!****************************************************************************************
678! 1) Flux calculation
679!****************************************************************************************
680! Suppose zero surface speed
681    u0(:)=0.0
682    v0(:)=0.0
683    u1_lay(:) = u1(:) - u0(:)
684    v1_lay(:) = v1(:) - v0(:)
685
686! set beta, cal, compute conduction fluxes inside ice/snow
687    slab_bilg(:)=0.
688    dif_grnd(:)=0.
689    beta(:) = 1.
690    DO i=1,knon
691    ki=knindex(i)
692        IF (snow(i).GT.snow_min) THEN
693            ! snow-layer heat capacity
694            cal(i)=2.*RCPD/(snow(i)*ice_cap)
695            ! snow conductive flux
696            f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
697            ! all shortwave flux absorbed
698            f_swpen=0.
699            ! bottom flux (ice conduction)
700            slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
701            ! update ice temperature
702            tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
703                     *(slab_bilg(ki)+f_cond)*dtime
704       ELSE ! bare ice
705            ! ice-layer heat capacity
706            cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
707            ! conductive flux
708            f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
709            ! penetrative shortwave flux...
710            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
711            slab_bilg(ki)=f_swpen-f_cond
712        END IF
713        radsol(i)=radsol(i)+f_cond-f_swpen
714    END DO
715    ! weight fluxes to ocean by sea ice fraction
716    slab_bilg(:)=slab_bilg(:)*fsic(:)
717
718! calcul_fluxs (sens, lat etc)
719    CALL calcul_fluxs(knon, is_sic, dtime, &
720        tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
721        precip_rain, precip_snow, snow, qsurf,  &
722        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
723        f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
724        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
725    DO i=1,knon
726        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
727    END DO
728
729! calcul_flux_wind
730    CALL calcul_flux_wind(knon, dtime, &
731         u0, v0, u1, v1, gustiness, cdragm, &
732         AcoefU, AcoefV, BcoefU, BcoefV, &
733         p1lay, temp_air, &
734         flux_u1, flux_v1)
735
736!****************************************************************************************
737! 2) Update snow and ice surface
738!****************************************************************************************
739! snow precip
740    DO i=1,knon
741        ki=knindex(i)
742        IF (precip_snow(i) > 0.) THEN
743            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
744        END IF
745! snow and ice sublimation
746        IF (evap(i) > 0.) THEN
747           snow_evap = MIN (snow(i) / dtime, evap(i))
748           snow(i) = snow(i) - snow_evap * dtime
749           snow(i) = MAX(0.0, snow(i))
750           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
751        ENDIF
752! Melt / Freeze snow from above if Tsurf>0
753        IF (tsurf_new(i).GT.t_melt) THEN
754            ! energy available for melting snow (in kg of melted snow /m2)
755            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
756               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
757            ! remove snow
758            IF (snow(i).GT.e_melt) THEN
759                snow(i)=snow(i)-e_melt
760                tsurf_new(i)=t_melt
761            ELSE ! all snow is melted
762                ! add remaining heat flux to ice
763                e_melt=e_melt-snow(i)
764                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
765                tsurf_new(i)=tice(ki)
766            END IF
767        END IF
768! melt ice from above if Tice>0
769        IF (tice(ki).GT.t_melt) THEN
770            ! quantity of ice melted (kg/m2)
771            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
772             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
773            ! melt from above, height only
774            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
775            e_melt=e_melt-dhsic
776            IF (e_melt.GT.0) THEN
777            ! lateral melt if ice too thin
778            dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
779            ! if all melted add remaining heat to ocean
780            e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
781            slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
782            ! update height and fraction
783            fsic(ki)=fsic(ki)-dfsic
784            END IF
785            seaice(ki)=seaice(ki)-dhsic
786            ! surface temperature at melting point
787            tice(ki)=t_melt
788            tsurf_new(i)=t_melt
789        END IF
790        ! convert snow to ice if below floating line
791        h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den
792        IF (h_test.GT.0.) THEN !snow under water
793            ! extra snow converted to ice (with added frozen sea water)
794            dhsic=h_test/(sea_den-ice_den+sno_den)
795            seaice(ki)=seaice(ki)+dhsic
796            snow(i)=snow(i)-dhsic*sno_den/ice_den
797            ! available energy (freeze sea water + bring to tice)
798            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
799                   ice_cap/2.*(t_freeze-tice(ki)))
800            ! update ice temperature
801            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
802        END IF
803    END DO
804
805! New albedo
806    DO i=1,knon
807        ki=knindex(i)
808       ! snow albedo: update snow age
809        IF (snow(i).GT.0.0001) THEN
810             agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
811                         * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.)
812        ELSE
813            agesno(i)=0.0
814        END IF
815        ! snow albedo
816        alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
817        ! ice albedo (varies with ice tkickness and temp)
818        alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
819        IF (tice(ki).GT.t_freeze-0.01) THEN
820            alb_ice=MIN(alb_ice,alb_ice_wet)
821        ELSE
822            alb_ice=MIN(alb_ice,alb_ice_dry)
823        END IF
824        ! pond albedo
825        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
826        ! pond fraction
827        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
828        ! snow fraction
829        frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
830        ! ice fraction
831        frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
832        ! total albedo
833        alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
834    END DO
835    alb2_new(:) = alb1_new(:)
836
837!****************************************************************************************
838! 3) Recalculate new ocean temperature (add fluxes below ice)
839!    Melt / freeze from below
840!***********************************************o*****************************************
841    !cumul fluxes
842    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
843    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
844        ! Add cumulated surface fluxes
845        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
846        DO i=1,knon
847            ki=knindex(i)
848            ! split lateral/top melt-freeze
849            frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
850            IF (tslab(ki,1).LE.t_freeze) THEN
851               ! ****** Form new ice from below *******
852               ! quantity of new ice
853                e_melt=(t_freeze-tslab(ki,1))/cyang &
854                       /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
855               ! first increase height to h_thin
856               dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
857               seaice(ki)=dhsic+seaice(ki)
858               e_melt=e_melt-fsic(ki)*dhsic
859               IF (e_melt.GT.0.) THEN
860               ! frac_mf fraction used for lateral increase
861               dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
862               fsic(ki)=fsic(ki)+dfsic
863               e_melt=e_melt-dfsic*seaice(ki)
864               ! rest used to increase height
865               seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
866               END IF
867               tslab(ki,1)=t_freeze
868           ELSE ! slab temperature above freezing
869               ! ****** melt ice from below *******
870               ! quantity of melted ice
871               e_melt=(tslab(ki,1)-t_freeze)/cyang &
872                       /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
873               ! first decrease height to h_thick
874               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
875               seaice(ki)=seaice(ki)-dhsic
876               e_melt=e_melt-fsic(ki)*dhsic
877               IF (e_melt.GT.0) THEN
878               ! frac_mf fraction used for height decrease
879               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
880               seaice(ki)=seaice(ki)-dhsic
881               e_melt=e_melt-fsic(ki)*dhsic
882               ! rest used to decrease fraction (up to 0!)
883               dfsic=MIN(fsic(ki),e_melt/seaice(ki))
884               ! keep remaining in ocean
885               e_melt=e_melt-dfsic*seaice(ki)
886               END IF
887               tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang
888               fsic(ki)=fsic(ki)-dfsic
889           END IF
890        END DO
891        bilg_cum(:)=0.
892    END IF ! coupling time
893   
894    !tests ice fraction
895    WHERE (fsic.LT.ice_frac_min)
896        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
897        tice=t_melt
898        fsic=0.
899        seaice=0.
900    END WHERE
901
902  END SUBROUTINE ocean_slab_ice
903!
904!****************************************************************************************
905!
906  SUBROUTINE ocean_slab_final
907
908!****************************************************************************************
909! Deallocate module variables
910!****************************************************************************************
911    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
912    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
913    IF (ALLOCATED(tice)) DEALLOCATE(tice)
914    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
915    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
916    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
917    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
918    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
919    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
920    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
921    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
922    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
923
924  END SUBROUTINE ocean_slab_final
925
926END MODULE ocean_slab_mod
Note: See TracBrowser for help on using the repository browser.