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

Last change on this file since 887 was 882, checked in by Laurent Fairhead, 18 years ago

Modifs pour intégrer le 1D 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: 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, qsurf, &
72       agesno, &
73       evap, fluxsens, fluxlat, &
74       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
75!
76! This subroutine treats the "open ocean", all grid points that are not entierly covered
77! by ice.
78! The routine reads data from climatologie file and does some calculations at the
79! surface.
80!
81    INCLUDE "indicesol.h"
82    INCLUDE "YOMCST.h"
83
84! Input arguments
85!****************************************************************************************
86    INTEGER, INTENT(IN)                      :: itime, jour, knon
87    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
88    LOGICAL, INTENT(IN)                      :: debut
89    REAL, INTENT(IN)                         :: dtime
90    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
91    REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
92    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
93    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
94    REAL, DIMENSION(klon), INTENT(IN)        :: petAcoef, peqAcoef
95    REAL, DIMENSION(klon), INTENT(IN)        :: petBcoef, peqBcoef
96    REAL, DIMENSION(klon), INTENT(IN)        :: ps
97    REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
98
99! In/Output arguments
100!****************************************************************************************
101    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
102    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
103    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
104 
105! Output arguments
106!****************************************************************************************
107    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
108    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
109    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
110    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
111    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
112
113! Local variables
114!****************************************************************************************
115    INTEGER                     :: i
116    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
117    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
118    LOGICAL                     :: check=.FALSE.
119    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
120
121!****************************************************************************************
122! Start calculation
123!****************************************************************************************
124    IF (check) WRITE(*,*)' Entering ocean_forced_noice'
125
126!****************************************************************************************
127! 1)   
128! Read from climatologie file SST and fraction of sub-surfaces
129!
130!****************************************************************************************
131! Get from file tsurf_lim and pctsrf_lim
132    CALL interfoce_lim(itime, dtime, jour, &
133         knon, knindex, &
134         debut, &
135         tsurf_lim, pctsrf_lim)
136   
137!****************************************************************************************
138! 2)
139! Flux calculation
140!
141!****************************************************************************************
142! Set some variables for calcul_fluxs
143    cal = 0.
144    beta = 1.
145    dif_grnd = 0.
146    alb_neig(:) = 0.
147    agesno(:) = 0.
148 
149! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
150    CALL calcul_fluxs(knon, is_oce, dtime, &
151         tsurf_lim, p1lay, cal, beta, tq_cdrag, ps, &
152         precip_rain, precip_snow, snow, qsurf,  &
153         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
154         petAcoef, peqAcoef, petBcoef, peqBcoef, &
155         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
156
157!****************************************************************************************
158! 3)
159! Calculate tmp_flux_o
160!
161!****************************************************************************************
162!IM: flux ocean-atmosphere utile pour le "slab" ocean
163! The flux are written to hist file
164    tmp_flux_o(:) = 0.0
165    DO i=1, knon
166       zx_sl(i) = RLVTT
167       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
168
169!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
170!       flux_o(i) = fluxsens(i) + fluxlat(i)
171       IF (pctsrf_lim(knindex(i),is_oce) .GT. epsfra) THEN
172          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
173       ENDIF
174    ENDDO
175
176!****************************************************************************************
177! 4)
178! Return the new values for the ocean fraction
179!
180!****************************************************************************************
181
182    pctsrf_oce(:) = pctsrf_lim(:,is_oce)
183
184  END SUBROUTINE ocean_forced_noice
185!
186!****************************************************************************************
187!
188  SUBROUTINE ocean_forced_ice(itime, dtime, jour, knon, knindex, &
189       debut, &
190       tsurf_in, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
191       petAcoef, peqAcoef, petBcoef, peqBcoef, &
192       ps, u1_lay, v1_lay, &
193       radsol, snow, qsurf, qsol, agesno, &
194       tsoil, alblw, evap, fluxsens, fluxlat, &
195       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
196!
197! This subroutine treats the ocean where there is ice.
198! The routine reads data from climatologie file and does flux calculations at the
199! surface.
200!   
201    INCLUDE "indicesol.h"
202    INCLUDE "dimsoil.h"
203    INCLUDE "YOMCST.h"
204    INCLUDE "clesphys.h"
205
206! Input arguments
207!****************************************************************************************
208    INTEGER, INTENT(IN)                  :: itime, jour, knon
209    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
210    LOGICAL, INTENT(IN)                  :: debut
211    REAL, INTENT(IN)                     :: dtime
212    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
213    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
214    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
215    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
216    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
217    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
218    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
219    REAL, DIMENSION(klon), INTENT(IN)    :: ps
220    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
221
222
223! In/Output arguments
224!****************************************************************************************
225    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
226    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
227    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
228    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
229
230
231! Output arguments
232!****************************************************************************************
233    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
234    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
235    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
236    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_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
289       CALL soil(dtime, is_sic, knon,snow, tsurf_tmp, tsoil,soilcap, soilflux)
290       cal(1:knon) = RCPD / soilcap(1:knon)
291       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
292       dif_grnd = 1.0 / tau_gl
293    ELSE
294       dif_grnd = 1.0 / tau_gl
295       cal = RCPD * calice
296       WHERE (snow > 0.0) cal = RCPD * calsno
297    ENDIF
298
299    beta = 1.0
300    CALL calcul_fluxs(knon, is_sic, dtime, &
301         tsurf_tmp, p1lay, cal, beta, tq_cdrag, ps, &
302         precip_rain, precip_snow, snow, qsurf,  &
303         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
304         petAcoef, peqAcoef, petBcoef, peqBcoef, &
305         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
306
307!****************************************************************************************
308! 3)
309! Calculations due to snow and runoff
310!
311!****************************************************************************************
312    CALL fonte_neige( knon, is_sic, knindex, dtime, &
313         tsurf_tmp, precip_rain, precip_snow, &
314         snow, qsol, tsurf_new, evap)
315   
316! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
317!
318    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 
319
320    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
321
322    alb_new(:) = 0.0
323    DO i=1, knon
324       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
325       alb_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
326    ENDDO
327!!      alb_new(1 : knon) = 0.6
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!!$    z0_new = 0.002
347!!$    z0_new = SQRT(z0_new**2+rugoro**2)
348    alblw(1:knon) = alb_new(1:knon)
349
350!****************************************************************************************
351! 5)
352! Return the new values for the seaice fraction
353!
354!****************************************************************************************
355
356    pctsrf_sic(:) = pctsrf_lim(:,is_sic)
357
358
359  END SUBROUTINE ocean_forced_ice
360!
361!****************************************************************************************
362!
363  SUBROUTINE ocean_forced_get_vars(flux_o, flux_g)
364! Get some variables from module oceanforced.
365! This subroutine returns variables to a external routine
366
367    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
368    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g
369
370! Initialize the output variables
371    flux_o(:) = tmp_flux_o(:)
372    flux_g(:) = tmp_flux_g(:)
373
374  END SUBROUTINE ocean_forced_get_vars
375!
376!****************************************************************************************
377!
378END MODULE ocean_forced_mod
379
380
381
382
383
384
Note: See TracBrowser for help on using the repository browser.