source: LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_forced_mod.F90 @ 5453

Last change on this file since 5453 was 4013, checked in by lguez, 3 years ago

Sync latest trunk changes to Ocean_skin

  • 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 Id
File size: 13.6 KB
Line 
1!
2! $Id: ocean_forced_mod.F90 4013 2021-11-19 15:58:59Z aborella $
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  IMPLICIT NONE
10
11CONTAINS
12!
13!****************************************************************************************
14!
15  SUBROUTINE ocean_forced_noice( &
16       itime, dtime, jour, knon, knindex, &
17       p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
18       temp_air, spechum, &
19       AcoefH, AcoefQ, BcoefH, BcoefQ, &
20       AcoefU, AcoefV, BcoefU, BcoefV, &
21       ps, u1, v1, gustiness, tsurf_in, &
22       radsol, snow, agesno, &
23       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
24       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
25!
26! This subroutine treats the "open ocean", all grid points that are not entierly covered
27! by ice.
28! The routine receives data from climatologie file limit.nc and does some calculations at the
29! surface.
30!
31    USE dimphy
32    USE calcul_fluxs_mod
33    USE limit_read_mod
34    USE mod_grid_phy_lmdz
35    USE indice_sol_mod
36    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
37    use config_ocean_skin_m, only: activate_ocean_skin
38
39    INCLUDE "YOMCST.h"
40    INCLUDE "clesphys.h"
41    INCLUDE "flux_arp.h"
42
43! Input arguments
44!****************************************************************************************
45    INTEGER, INTENT(IN)                      :: itime, jour, knon
46    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
47    REAL, INTENT(IN)                         :: dtime
48    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
49    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh, cdragq, cdragm
50    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
51    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
52    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
53    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
54    REAL, DIMENSION(klon), INTENT(IN)        :: ps
55    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
56    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
57    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
58
59! In/Output arguments
60!****************************************************************************************
61    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
62    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
63    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
64 
65! Output arguments
66!****************************************************************************************
67    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
68    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
69    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
70    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
71    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
72    REAL, intent(out):: sens_prec_liq(:) ! (knon)
73
74! Local variables
75!****************************************************************************************
76    INTEGER                     :: i, j
77    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
78    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
79    REAL, DIMENSION(klon)       :: u0, v0
80    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
81    LOGICAL                     :: check=.FALSE.
82    REAL sens_prec_sol(knon)
83    REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol   
84
85!****************************************************************************************
86! Start calculation
87!****************************************************************************************
88    IF (check) WRITE(*,*)' Entering ocean_forced_noice'
89   
90!****************************************************************************************
91! 1)   
92! Read sea-surface temperature from file limit.nc
93!
94!****************************************************************************************
95!--sb:
96!!jyg    if (knon.eq.1) then ! single-column model
97    if (klon_glo.eq.1) then ! single-column model
98      ! EV: now surface Tin flux_arp.h
99      !CALL read_tsurf1d(knon,tsurf_lim) ! new
100       DO i = 1, knon
101        tsurf_lim(i) = tg
102       ENDDO
103
104    else ! GCM
105      CALL limit_read_sst(knon,knindex,tsurf_lim)
106    endif ! knon
107!sb--
108
109!****************************************************************************************
110! 2)
111! Flux calculation
112!
113!****************************************************************************************
114! Set some variables for calcul_fluxs
115    !cal = 0.
116    !beta = 1.
117    !dif_grnd = 0.
118   
119   
120    ! EV: use calbeta to calculate beta
121    ! Need to initialize qsurf for calbeta but it is not modified by this routine
122    qsurf(:)=0.
123    CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd)
124
125
126    alb_neig(:) = 0.
127    agesno(:) = 0.
128    lat_prec_liq = 0.; lat_prec_sol = 0.
129
130! Suppose zero surface speed
131    u0(:)=0.0
132    v0(:)=0.0
133    u1_lay(:) = u1(:) - u0(:)
134    v1_lay(:) = v1(:) - v0(:)
135
136! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
137    CALL calcul_fluxs(knon, is_oce, dtime, &
138         merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), p1lay, cal, &
139         beta, cdragh, cdragq, ps, &
140         precip_rain, precip_snow, snow, qsurf,  &
141         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
142         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
143         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
144         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
145    if (activate_ocean_skin == 2) tsurf_new = tsurf_lim
146
147    do j = 1, knon
148      i = knindex(j)
149      sens_prec_liq_o(i,1) = sens_prec_liq(j)
150      sens_prec_sol_o(i,1) = sens_prec_sol(j)
151      lat_prec_liq_o(i,1) = lat_prec_liq(j)
152      lat_prec_sol_o(i,1) = lat_prec_sol(j)
153    enddo
154
155
156! - Flux calculation at first modele level for U and V
157    CALL calcul_flux_wind(knon, dtime, &
158         u0, v0, u1, v1, gustiness, cdragm, &
159         AcoefU, AcoefV, BcoefU, BcoefV, &
160         p1lay, temp_air, &
161         flux_u1, flux_v1) 
162
163  END SUBROUTINE ocean_forced_noice
164!
165!***************************************************************************************
166!
167  SUBROUTINE ocean_forced_ice( &
168       itime, dtime, jour, knon, knindex, &
169       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
170       AcoefH, AcoefQ, BcoefH, BcoefQ, &
171       AcoefU, AcoefV, BcoefU, BcoefV, &
172       ps, u1, v1, gustiness, &
173       radsol, snow, qsol, agesno, tsoil, &
174       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
175       tsurf_new, dflux_s, dflux_l, rhoa)
176!
177! This subroutine treats the ocean where there is ice.
178! The routine reads data from climatologie file and does flux calculations at the
179! surface.
180!
181    USE dimphy
182    USE geometry_mod, ONLY: longitude,latitude
183    USE calcul_fluxs_mod
184    USE surface_data,     ONLY : calice, calsno
185    USE limit_read_mod
186    USE fonte_neige_mod,  ONLY : fonte_neige
187    USE indice_sol_mod
188    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
189
190!   INCLUDE "indicesol.h"
191    INCLUDE "dimsoil.h"
192    INCLUDE "YOMCST.h"
193    INCLUDE "clesphys.h"
194    INCLUDE "flux_arp.h"
195
196! Input arguments
197!****************************************************************************************
198    INTEGER, INTENT(IN)                  :: itime, jour, knon
199    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
200    REAL, INTENT(IN)                     :: dtime
201    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
202    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
203    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
204    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
205    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
206    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
207    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
208    REAL, DIMENSION(klon), INTENT(IN)    :: ps
209    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
210    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
211
212! In/Output arguments
213!****************************************************************************************
214    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
215    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
216    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
217    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
218
219! Output arguments
220!****************************************************************************************
221    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
222    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
223    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
224    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
225    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
226    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
227    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
228
229! Local variables
230!****************************************************************************************
231    LOGICAL                     :: check=.FALSE.
232    INTEGER                     :: i, j
233    REAL                        :: zfra
234    REAL, PARAMETER             :: t_grnd=271.35
235    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
236    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
237    REAL, DIMENSION(klon)       :: soilcap, soilflux
238    REAL, DIMENSION(klon)       :: u0, v0
239    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
240    REAL sens_prec_liq(knon), sens_prec_sol (knon)   
241    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol   
242
243
244!****************************************************************************************
245! Start calculation
246!****************************************************************************************
247    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon
248
249!****************************************************************************************
250! 1)
251! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
252!                    dflux_s, dflux_l and qsurf
253!****************************************************************************************
254
255    tsurf_tmp(:) = tsurf_in(:)
256
257! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal
258    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
259
260   
261    IF (soil_model) THEN
262! update tsoil and calculate soilcap and soilflux
263       CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, &
264        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil,soilcap, soilflux)
265       cal(1:knon) = RCPD / soilcap(1:knon)
266       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
267       dif_grnd = 1.0 / tau_gl
268    ELSE
269       dif_grnd = 1.0 / tau_gl
270       cal = RCPD * calice
271       WHERE (snow > 0.0) cal = RCPD * calsno
272    ENDIF
273
274!    beta = 1.0
275    lat_prec_liq = 0.; lat_prec_sol = 0.
276
277! Suppose zero surface speed
278    u0(:)=0.0
279    v0(:)=0.0
280    u1_lay(:) = u1(:) - u0(:)
281    v1_lay(:) = v1(:) - v0(:)
282    CALL calcul_fluxs(knon, is_sic, dtime, &
283         tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, &
284         precip_rain, precip_snow, snow, qsurf,  &
285         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
286         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
287         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
288         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
289    do j = 1, knon
290      i = knindex(j)
291      sens_prec_liq_o(i,2) = sens_prec_liq(j)
292      sens_prec_sol_o(i,2) = sens_prec_sol(j)
293      lat_prec_liq_o(i,2) = lat_prec_liq(j)
294      lat_prec_sol_o(i,2) = lat_prec_sol(j)
295    enddo
296
297! - Flux calculation at first modele level for U and V
298    CALL calcul_flux_wind(knon, dtime, &
299         u0, v0, u1, v1, gustiness, cdragm, &
300         AcoefU, AcoefV, BcoefU, BcoefV, &
301         p1lay, temp_air, &
302         flux_u1, flux_v1) 
303
304!****************************************************************************************
305! 2)
306! Calculations due to snow and runoff
307!
308!****************************************************************************************
309    CALL fonte_neige( knon, is_sic, knindex, dtime, &
310         tsurf_tmp, precip_rain, precip_snow, &
311         snow, qsol, tsurf_new, evap)
312   
313! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
314!
315    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 
316
317    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
318
319    alb1_new(:) = 0.0
320    DO i=1, knon
321       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
322       alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
323    ENDDO
324
325    alb2_new(:) = alb1_new(:)
326
327  END SUBROUTINE ocean_forced_ice
328
329!************************************************************************
330! 1D case
331!************************************************************************
332!  SUBROUTINE read_tsurf1d(knon,sst_out)
333!
334! This subroutine specifies the surface temperature to be used in 1D simulations
335!
336!      USE dimphy, ONLY : klon
337!
338!      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
339!      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
340!
341!       INTEGER :: i
342! COMMON defined in lmdz1d.F:
343!       real ts_cur
344!       common /sst_forcing/ts_cur
345!
346!       DO i = 1, knon
347!        sst_out(i) = ts_cur
348!       ENDDO
349!
350!      END SUBROUTINE read_tsurf1d
351!
352!
353!************************************************************************
354END MODULE ocean_forced_mod
355
356
357
358
359
360
Note: See TracBrowser for help on using the repository browser.