source: LMDZ4/trunk/libf/phytherm/ocean_forced_mod.F90 @ 1068

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

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