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

Last change on this file since 2195 was 2057, checked in by Ehouarn Millour, 10 years ago

Preparatory stuff to fix and improve the slab ocean model.
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: 12.3 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
10
[781]11  IMPLICIT NONE
[996]12  PRIVATE
[2057]13  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
[781]14
[2057]15  INTEGER, PRIVATE, SAVE                           :: cpl_pas
16  !$OMP THREADPRIVATE(cpl_pas)
17  REAL, PRIVATE, SAVE                              :: cyang
18  !$OMP THREADPRIVATE(cyang)
19  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
20  !$OMP THREADPRIVATE(slabh)
21  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
22  !$OMP THREADPRIVATE(tslab)
23  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
24  !$OMP THREADPRIVATE(pctsrf)
25  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
26  !$OMP THREADPRIVATE(slab_bils)
27  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
28  !$OMP THREADPRIVATE(bils_cum)
29
[781]30CONTAINS
31!
32!****************************************************************************************
33!
[2057]34  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
35  !, seaice_rst etc
[781]36
[2057]37    use IOIPSL
38
39    INCLUDE "iniprint.h"
40    ! For ok_xxx vars (Ekman...)
41    INCLUDE "clesphys.h"
42
43    ! Input variables
44!****************************************************************************************
45    REAL, INTENT(IN)                         :: dtime
46! Variables read from restart file
47    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
48
49! Local variables
50!****************************************************************************************
51    INTEGER                :: error
52    CHARACTER (len = 80)   :: abort_message
53    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
54
55!****************************************************************************************
56! Allocate surface fraction read from restart file
57!****************************************************************************************
58    ALLOCATE(pctsrf(klon,nbsrf), stat = error)
59    IF (error /= 0) THEN
60       abort_message='Pb allocation tmp_pctsrf_slab'
61       CALL abort_gcm(modname,abort_message,1)
62    ENDIF
63    pctsrf(:,:) = pctsrf_rst(:,:)
64
65!****************************************************************************************
66! Allocate local variables
67!****************************************************************************************
68    ALLOCATE(slab_bils(klon), stat = error)
69    IF (error /= 0) THEN
70       abort_message='Pb allocation slab_bils'
71       CALL abort_gcm(modname,abort_message,1)
72    ENDIF
73    slab_bils(:) = 0.0   
74    ALLOCATE(bils_cum(klon), stat = error)
75    IF (error /= 0) THEN
76       abort_message='Pb allocation slab_bils_cum'
77       CALL abort_gcm(modname,abort_message,1)
78    ENDIF
79    bils_cum(:) = 0.0   
80
81! Layer thickness
82    ALLOCATE(slabh(nslay), stat = error)
83    IF (error /= 0) THEN
84       abort_message='Pb allocation slabh'
85       CALL abort_gcm(modname,abort_message,1)
86    ENDIF
87    slabh(1)=50.
88! cyang = 1/heat capacity of top layer (rho.c.H)
89    cyang=1/(slabh(1)*4.228e+06)
90
91! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
92! pour un calcul à chaque pas de temps, cpl_pas=1
93    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
94    CALL getin('cpl_pas',cpl_pas)
95    print *,'cpl_pas',cpl_pas
96 END SUBROUTINE ocean_slab_init
97!
98!****************************************************************************************
99!
100  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
101
[996]102    USE limit_read_mod
[1067]103    USE surface_data
[1785]104
[996]105!    INCLUDE "clesphys.h"
[781]106
[996]107! Arguments
[781]108!****************************************************************************************
[996]109    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
110    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
111    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
[2057]112    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
[996]113    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
[781]114
115! Local variables
116!****************************************************************************************
117
[2057]118    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
119       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
[996]120    ELSE
[2057]121       pctsrf_chg(:,:)=pctsrf(:,:)
122       is_modified=.TRUE.
[996]123    END IF
[781]124
[996]125  END SUBROUTINE ocean_slab_frac
[781]126!
127!****************************************************************************************
128!
129  SUBROUTINE ocean_slab_noice( &
[996]130       itime, dtime, jour, knon, knindex, &
[1067]131       p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
132       AcoefH, AcoefQ, BcoefH, BcoefQ, &
133       AcoefU, AcoefV, BcoefU, BcoefV, &
134       ps, u1, v1, tsurf_in, &
[888]135       radsol, snow, agesno, &
[1067]136       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[2057]137       tsurf_new, dflux_s, dflux_l, qflux)
[1067]138   
139    USE calcul_fluxs_mod
[2057]140    USE surface_data
[1785]141
[781]142    INCLUDE "iniprint.h"
143
144! Input arguments
145!****************************************************************************************
[996]146    INTEGER, INTENT(IN)                  :: itime
147    INTEGER, INTENT(IN)                  :: jour
[781]148    INTEGER, INTENT(IN)                  :: knon
149    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
150    REAL, INTENT(IN)                     :: dtime
151    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
[1067]152    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
[781]153    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
154    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
[1067]155    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
156    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
[781]157    REAL, DIMENSION(klon), INTENT(IN)    :: ps
[1067]158    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
[996]159    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
[781]160
161! In/Output arguments
162!****************************************************************************************
163    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
164    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
165    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
166   
167! Output arguments
168!****************************************************************************************
169    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
170    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
[1067]171    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
[781]172    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
173    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
[2057]174    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
[781]175
176! Local variables
177!****************************************************************************************
[2057]178    INTEGER               :: i,ki
[996]179    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
[2057]180    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
[1067]181    REAL, DIMENSION(klon) :: u0, v0
182    REAL, DIMENSION(klon) :: u1_lay, v1_lay
[781]183
184!****************************************************************************************
[996]185! 1) Flux calculation
186!
187!****************************************************************************************
188    cal(:)      = 0.
189    beta(:)     = 1.
190    dif_grnd(:) = 0.
191    agesno(:)   = 0.
[781]192   
[1067]193! Suppose zero surface speed
194    u0(:)=0.0
195    v0(:)=0.0
196    u1_lay(:) = u1(:) - u0(:)
197    v1_lay(:) = v1(:) - v0(:)
198
[781]199    CALL calcul_fluxs(knon, is_oce, dtime, &
[1067]200         tsurf_in, p1lay, cal, beta, cdragh, ps, &
[781]201         precip_rain, precip_snow, snow, qsurf,  &
202         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
[1067]203         AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]204         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
205
[1067]206! - Flux calculation at first modele level for U and V
207    CALL calcul_flux_wind(knon, dtime, &
208         u0, v0, u1, v1, cdragm, &
209         AcoefU, AcoefV, BcoefU, BcoefV, &
210         p1lay, temp_air, &
211         flux_u1, flux_v1) 
212
[2057]213! Accumulate total fluxes locally
214    slab_bils(:)=0.
215    DO i=1,knon
216        ki=knindex(i)
217        slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
218        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
219! Also taux, tauy, saved vars...
220    END DO
221
[781]222!****************************************************************************************
[2057]223! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
[781]224!
225!****************************************************************************************
[2057]226    lmt_bils(:)=0.
227    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst)  ! global pour un processus
228    ! lmt_bils and diff_sst saved by limit_slab
229    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
230    ! qflux = total QFlux correction (in W/m2)
[781]231
232!****************************************************************************************
[996]233! 3) Recalculate new temperature
[781]234!
[2057]235!***********************************************o*****************************************
236    tsurf_new=tsurf_in
237    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
238        ! Compute transport
239        ! Add read QFlux and SST tendency
240        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
241        ! Add cumulated surface fluxes
242        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
243        ! Update surface temperature
244        SELECT CASE(version_ocean)
245        CASE('sicNO')
246            DO i=1,knon
247                ki=knindex(i)
248                tsurf_new(i)=tslab(ki,1)
249            END DO
250        CASE('sicOBS') ! check for sea ice or tsurf below freezing
251            DO i=1,knon
252                ki=knindex(i)
253                IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
254                    tsurf_new(i)=t_freeze
255                    tslab(ki,1)=t_freeze
256                ELSE
257                    tsurf_new(i)=tslab(ki,1)
258                END IF
259            END DO
260        CASE('sicINT')
261            DO i=1,knon
262                ki=knindex(i)
263                IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
264                    IF (tslab(ki,1).GT.t_freeze) THEN
265                        tsurf_new(i)=tslab(ki,1)
266                    ELSE
267                        tsurf_new(i)=t_freeze
268                        ! Call new ice routine
269                        tslab(ki,1)=t_freeze
270                    END IF
271                ELSE ! ice present, tslab update completed in slab_ice
272                    tsurf_new(i)=t_freeze
273                END IF !ice free
274            END DO
275        END SELECT
276        bils_cum(:)=0.0! clear cumulated fluxes
277    END IF ! coupling time
278  END SUBROUTINE ocean_slab_noice
279!
[781]280!****************************************************************************************
[2057]281!
282!  SUBROUTINE ocean_slab_ice(   &
283!       itime, dtime, jour, knon, knindex, &
284!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
285!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
286!       AcoefU, AcoefV, BcoefU, BcoefV, &
287!       ps, u1, v1, &
288!       radsol, snow, qsurf, qsol, agesno, tsoil, &
289!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
290!       tsurf_new, dflux_s, dflux_l)
291!
292!****************************************************************************************
293! 1) Flux calculation
294!****************************************************************************************
295! set beta, cal etc. depends snow / ice surf ?
296! calcul_fluxs (sens, lat etc)
297! calcul_flux_wind
[781]298
[2057]299!****************************************************************************************
300! 2) Update surface
301!****************************************************************************************
302! neige, fonte
303! flux glace-ocean
304! update temperature
305! neige precip, evap
306! Melt snow & ice from above
307! New albedo
308
309!****************************************************************************************
310! 3) Recalculate new ocean temperature
311!    Melt / freeze from below
312!***********************************************o*****************************************
313
314
315!  END SUBROUTINE ocean_slab_ice
[781]316!
317!****************************************************************************************
318!
[2057]319  SUBROUTINE ocean_slab_final
320  !, seaice_rst etc
321
322    ! For ok_xxx vars (Ekman...)
323    INCLUDE "clesphys.h"
324
325!****************************************************************************************
326! Deallocate module variables
327!
328!****************************************************************************************
329    IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
330    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
331
332  END SUBROUTINE ocean_slab_final
333
[781]334END MODULE ocean_slab_mod
Note: See TracBrowser for help on using the repository browser.