source: LMDZ6/branches/LMDZ-tracers/libf/phylmd/ocean_slab_mod.F90 @ 3871

Last change on this file since 3871 was 3851, checked in by dcugnet, 4 years ago

Update the branch to the current trunk.

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