source: LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90 @ 922

Last change on this file since 922 was 888, checked in by Laurent Fairhead, 17 years ago

Modifications sur l'albedo JG
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 14.0 KB
Line 
1!
2! $Header$
3!
4MODULE ocean_forced_mod
5!
6! This module is used for both the sub-surfaces ocean and sea-ice for the case of a
7! forced ocean,  "ocean=force".
8!
9  USE surface_data,     ONLY : calice, calsno, tau_gl
10  USE fonte_neige_mod,  ONLY : fonte_neige
11  USE calcul_fluxs_mod, ONLY : calcul_fluxs
12  USE dimphy
13
14  IMPLICIT NONE
15
16  REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE   :: tmp_flux_o, tmp_flux_g
17  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
18
19CONTAINS
20!
21!****************************************************************************************
22!
23  SUBROUTINE ocean_forced_init
24! Allocate fields needed for this module
25!   
26    INTEGER              :: error
27    CHARACTER (len = 80) :: abort_message
28    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
29!****************************************************************************************
30
31    ALLOCATE(tmp_flux_o(1:klon), stat = error)
32    IF (error /= 0) THEN
33       abort_message='Pb allocation tmp_flux_o'
34       CALL abort_gcm(modname,abort_message,1)
35    ENDIF
36
37    ALLOCATE(tmp_flux_g(1:klon), stat = error)
38    IF (error /= 0) THEN
39       abort_message='Pb allocation tmp_flux_g'
40       CALL abort_gcm(modname,abort_message,1)
41    ENDIF   
42
43
44  END SUBROUTINE ocean_forced_init
45!
46!****************************************************************************************
47!****************************************************************************************
48!
49  SUBROUTINE ocean_forced_final
50! Allocate fields needed for this module
51!   
52    INTEGER              :: error
53    CHARACTER (len = 80) :: abort_message
54    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
55!****************************************************************************************
56
57    DEALLOCATE(tmp_flux_o)
58    DEALLOCATE(tmp_flux_g)
59
60
61  END SUBROUTINE ocean_forced_final
62!
63!****************************************************************************************
64!
65  SUBROUTINE ocean_forced_noice(itime, dtime, jour, knon, knindex, &
66       debut, &
67       p1lay, tq_cdrag, precip_rain, precip_snow, &
68       temp_air, spechum, &
69       petAcoef, peqAcoef, petBcoef, peqBcoef, &
70       ps, u1_lay, v1_lay, &
71       radsol, snow, agesno, &
72       qsurf, evap, fluxsens, fluxlat, &
73       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
74!
75! This subroutine treats the "open ocean", all grid points that are not entierly covered
76! by ice.
77! The routine reads data from climatologie file and does some calculations at the
78! surface.
79!
80    INCLUDE "indicesol.h"
81    INCLUDE "YOMCST.h"
82
83! Input arguments
84!****************************************************************************************
85    INTEGER, INTENT(IN)                      :: itime, jour, knon
86    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
87    LOGICAL, INTENT(IN)                      :: debut
88    REAL, INTENT(IN)                         :: dtime
89    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
90    REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
91    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
92    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
93    REAL, DIMENSION(klon), INTENT(IN)        :: petAcoef, peqAcoef
94    REAL, DIMENSION(klon), INTENT(IN)        :: petBcoef, peqBcoef
95    REAL, DIMENSION(klon), INTENT(IN)        :: ps
96    REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
97
98! In/Output arguments
99!****************************************************************************************
100    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
101    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
102    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
103 
104! Output arguments
105!****************************************************************************************
106    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
107    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
108    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
109    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
110    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
111
112! Local variables
113!****************************************************************************************
114    INTEGER                     :: i
115    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
116    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
117    LOGICAL                     :: check=.FALSE.
118    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
119
120!****************************************************************************************
121! Start calculation
122!****************************************************************************************
123    IF (check) WRITE(*,*)' Entering ocean_forced_noice'
124
125!****************************************************************************************
126! 1)   
127! Read from climatologie file SST and fraction of sub-surfaces
128!
129!****************************************************************************************
130! Get from file tsurf_lim and pctsrf_lim
131    CALL interfoce_lim(itime, dtime, jour, &
132         knon, knindex, &
133         debut, &
134         tsurf_lim, pctsrf_lim)
135   
136!****************************************************************************************
137! 2)
138! Flux calculation
139!
140!****************************************************************************************
141! Set some variables for calcul_fluxs
142    cal = 0.
143    beta = 1.
144    dif_grnd = 0.
145    alb_neig(:) = 0.
146    agesno(:) = 0.
147 
148! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
149    CALL calcul_fluxs(knon, is_oce, dtime, &
150         tsurf_lim, p1lay, cal, beta, tq_cdrag, ps, &
151         precip_rain, precip_snow, snow, qsurf,  &
152         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
153         petAcoef, peqAcoef, petBcoef, peqBcoef, &
154         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
155
156!****************************************************************************************
157! 3)
158! Calculate tmp_flux_o
159!
160!****************************************************************************************
161!IM: flux ocean-atmosphere utile pour le "slab" ocean
162! The flux are written to hist file
163    tmp_flux_o(:) = 0.0
164    DO i=1, knon
165       zx_sl(i) = RLVTT
166       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
167
168!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
169!       flux_o(i) = fluxsens(i) + fluxlat(i)
170       IF (pctsrf_lim(knindex(i),is_oce) .GT. epsfra) THEN
171          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
172       ENDIF
173    ENDDO
174
175!****************************************************************************************
176! 4)
177! Return the new values for the ocean fraction
178!
179!****************************************************************************************
180
181    pctsrf_oce(:) = pctsrf_lim(:,is_oce)
182
183  END SUBROUTINE ocean_forced_noice
184!
185!****************************************************************************************
186!
187  SUBROUTINE ocean_forced_ice(itime, dtime, jour, knon, knindex, &
188       debut, &
189       tsurf_in, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
190       petAcoef, peqAcoef, petBcoef, peqBcoef, &
191       ps, u1_lay, v1_lay, &
192       radsol, snow, qsol, agesno, tsoil, &
193       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
194       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
195!
196! This subroutine treats the ocean where there is ice.
197! The routine reads data from climatologie file and does flux calculations at the
198! surface.
199!   
200    INCLUDE "indicesol.h"
201    INCLUDE "dimsoil.h"
202    INCLUDE "YOMCST.h"
203    INCLUDE "clesphys.h"
204
205! Input arguments
206!****************************************************************************************
207    INTEGER, INTENT(IN)                  :: itime, jour, knon
208    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
209    LOGICAL, INTENT(IN)                  :: debut
210    REAL, INTENT(IN)                     :: dtime
211    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
212    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
213    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
214    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
215    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
216    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
217    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
218    REAL, DIMENSION(klon), INTENT(IN)    :: ps
219    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
220
221
222! In/Output arguments
223!****************************************************************************************
224    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
225    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
226    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
227    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
228
229
230! Output arguments
231!****************************************************************************************
232    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
233    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
234    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
235    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
236    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
237    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
238    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
239
240
241! Local variables
242!****************************************************************************************
243    LOGICAL                     :: check=.FALSE.
244    INTEGER                     :: i
245    REAL                        :: zfra
246    REAL, PARAMETER             :: t_grnd=271.35
247    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
248    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
249    REAL, DIMENSION(klon)       :: soilcap, soilflux
250    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
251
252!****************************************************************************************
253! Start calculation
254!****************************************************************************************
255    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon
256
257!****************************************************************************************
258! 1)   
259! Read from climatologie file SST and fraction of sub-surfaces
260!
261!****************************************************************************************
262    CALL interfoce_lim(itime, dtime, jour, &
263         knon, knindex, &
264         debut, &
265         tsurf_tmp, pctsrf_lim)
266 
267    DO i = 1, knon
268       tsurf_tmp(i) = tsurf_in(i)
269       IF (pctsrf_lim(knindex(i),is_sic) < EPSFRA) THEN
270          snow(i) = 0.0
271          tsurf_tmp(i) = RTT - 1.8
272          IF (soil_model) tsoil(i,:) = RTT -1.8
273       ENDIF
274    ENDDO
275
276!****************************************************************************************
277! 2)
278! Flux calculation : tsurf_new, evap, fluxlat, fluxsens,
279!                    dflux_s, dflux_l and qsurf
280!****************************************************************************************
281
282! calculate the parameters cal, beta, capsol and dif_grnd
283
284    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
285   
286    IF (soil_model) THEN
287! update tsoil and calculate soilcap and soilflux
288       CALL soil(dtime, is_sic, knon,snow, tsurf_tmp, tsoil,soilcap, soilflux)
289       cal(1:knon) = RCPD / soilcap(1:knon)
290       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
291       dif_grnd = 1.0 / tau_gl
292    ELSE
293       dif_grnd = 1.0 / tau_gl
294       cal = RCPD * calice
295       WHERE (snow > 0.0) cal = RCPD * calsno
296    ENDIF
297
298    beta = 1.0
299    CALL calcul_fluxs(knon, is_sic, dtime, &
300         tsurf_tmp, p1lay, cal, beta, tq_cdrag, ps, &
301         precip_rain, precip_snow, snow, qsurf,  &
302         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
303         petAcoef, peqAcoef, petBcoef, peqBcoef, &
304         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
305
306!****************************************************************************************
307! 3)
308! Calculations due to snow and runoff
309!
310!****************************************************************************************
311    CALL fonte_neige( knon, is_sic, knindex, dtime, &
312         tsurf_tmp, precip_rain, precip_snow, &
313         snow, qsol, tsurf_new, evap)
314   
315! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
316!
317    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 
318
319    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
320
321    alb1_new(:) = 0.0
322    DO i=1, knon
323       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
324       alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
325    ENDDO
326
327    alb2_new(:) = alb1_new(:)
328
329!****************************************************************************************
330! 4)
331! Calculate tmp_flux_g
332!
333!****************************************************************************************
334!
335    tmp_flux_g(:) = 0.0
336    DO i = 1, knon
337!IM: faire dependre le coefficient de conduction de la glace de mer
338!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
339!    actuel correspond a 3m de glace de mer, cf. L.Li
340!
341       IF (cal(i) .GT. 1.0e-15 .AND. pctsrf_lim(knindex(i),is_sic) .GT. epsfra) &
342            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * dif_grnd(i) *RCPD/cal(i)
343       
344    ENDDO
345
346
347!****************************************************************************************
348! 5)
349! Return the new values for the seaice fraction
350!
351!****************************************************************************************
352
353    pctsrf_sic(:) = pctsrf_lim(:,is_sic)
354
355  END SUBROUTINE ocean_forced_ice
356!
357!****************************************************************************************
358!
359  SUBROUTINE ocean_forced_get_vars(flux_o, flux_g)
360! Get some variables from module oceanforced.
361! This subroutine returns variables to a external routine
362
363    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
364    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g
365
366! Initialize the output variables
367    flux_o(:) = tmp_flux_o(:)
368    flux_g(:) = tmp_flux_g(:)
369
370  END SUBROUTINE ocean_forced_get_vars
371!
372!****************************************************************************************
373!
374END MODULE ocean_forced_mod
375
376
377
378
379
380
Note: See TracBrowser for help on using the repository browser.