source: LMDZ6/trunk/libf/phylmdiso/ocean_forced_mod.F90 @ 3929

Last change on this file since 3929 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 18.5 KB
RevLine 
[3927]1!
2! $Id: ocean_forced_mod.F90 3327 2018-05-16 13:58:59Z musat $
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, &
22       radsol, snow, agesno, &
23       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
24       tsurf_new, dflux_s, dflux_l &
25#ifdef ISO
26            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
27            xtsnow,xtevap,h1 & 
28#endif           
29            )
30!
31! This subroutine treats the "open ocean", all grid points that are not entierly covered
32! by ice.
33! The routine receives data from climatologie file limit.nc and does some calculations at the
34! surface.
35!
36    USE dimphy
37    USE calcul_fluxs_mod
38    USE limit_read_mod
39    USE mod_grid_phy_lmdz
40    USE indice_sol_mod
41    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
42#ifdef ISO
43  USE infotrac_phy, ONLY: ntraciso,niso
44    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, &
45&       calcul_iso_surf_sic_vectall   
46#ifdef ISOVERIF
47    USE isotopes_mod, ONLY: iso_eau,ridicule
48    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
49    USE isotopes_verif_mod
50#endif
51#endif
52
53    INCLUDE "YOMCST.h"
54    INCLUDE "clesphys.h"
55
56
57! Input arguments
58!****************************************************************************************
59    INTEGER, INTENT(IN)                      :: itime, jour, knon
60    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
61    REAL, INTENT(IN)                         :: dtime
62    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
63    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh, cdragq, cdragm
64    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
65    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
66    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
67    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
68    REAL, DIMENSION(klon), INTENT(IN)        :: ps
69    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
70#ifdef ISO
71    REAL, DIMENSION(ntraciso,klon), INTENT(IN)    :: xtprecip_rain, xtprecip_snow
72    REAL, DIMENSION(ntraciso,klon), INTENT(IN)    :: xtspechum
73    real, dimension(klon), intent(IN) :: rlat
74#endif
75
76! In/Output arguments
77!****************************************************************************************
78    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
79    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
80    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
81#ifdef ISO     
82    REAL, DIMENSION(niso,klon), INTENT(IN)    :: xtsnow
83    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: Roce
84#endif
85
86! Output arguments
87!****************************************************************************************
88    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
89    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
90    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
91    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
92    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
93#ifdef ISO     
94    REAL, DIMENSION(ntraciso,klon), INTENT(OUT)    :: xtevap ! isotopes in evaporation flux
95    REAL, DIMENSION(klon), INTENT(out)    :: h1 ! just a diagnostic, not useful for the simulation
96#endif
97
98! Local variables
99!****************************************************************************************
100    INTEGER                     :: i, j
101    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
102    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
103    REAL, DIMENSION(klon)       :: u0, v0
104    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
105    LOGICAL                     :: check=.FALSE.
106    REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol   
107    REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol   
108
109#ifdef ISO   
110      real, parameter :: t_coup = 273.15     
111#endif
112
113!****************************************************************************************
114! Start calculation
115!****************************************************************************************
116    IF (check) WRITE(*,*)' Entering ocean_forced_noice'
117   
118#ifdef ISO
119#ifdef ISOVERIF
120        do i=1,knon
121          if (iso_eau.gt.0) then         
122              call iso_verif_egalite_choix(xtspechum(iso_eau,i), &
123     &                  spechum(i),'ocean_forced_mod 111', &
124     &                  errmax,errmaxrel)     
125              call iso_verif_egalite_choix(snow(i), &
126     &          xtsnow(iso_eau,i),'ocean_forced_mod 117', &
127     &          errmax,errmaxrel)
128           endif !if (iso_eau.gt.0) then
129        enddo !do i=1,knon
130#endif     
131#endif
132
133!****************************************************************************************
134! 1)   
135! Read sea-surface temperature from file limit.nc
136!
137!****************************************************************************************
138!--sb:
139!!jyg    if (knon.eq.1) then ! single-column model
140    if (klon_glo.eq.1) then ! single-column model
141      CALL read_tsurf1d(knon,tsurf_lim) ! new
142#ifdef ISO
143      write(*,*) 'ocean_forced_mod 143: isotopes pas prévus ici'
144      stop
145#endif
146    else ! GCM
147      CALL limit_read_sst(knon,knindex,tsurf_lim &
148#ifdef ISO
149     &  ,Roce,rlat   &
150#endif     
151     &   )
152    endif ! knon
153!sb--
154
155!****************************************************************************************
156! 2)
157! Flux calculation
158!
159!****************************************************************************************
160! Set some variables for calcul_fluxs
161    cal = 0.
162    beta = 1.
163    dif_grnd = 0.
164    alb_neig(:) = 0.
165    agesno(:) = 0.
166    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
167
168! Suppose zero surface speed
169    u0(:)=0.0
170    v0(:)=0.0
171    u1_lay(:) = u1(:) - u0(:)
172    v1_lay(:) = v1(:) - v0(:)
173
174! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
175    CALL calcul_fluxs(knon, is_oce, dtime, &
176         tsurf_lim, p1lay, cal, beta, cdragh, cdragq, ps, &
177         precip_rain, precip_snow, snow, qsurf,  &
178         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
179         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
180         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
181         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
182
183    do j = 1, knon
184      i = knindex(j)
185      sens_prec_liq_o(i,1) = sens_prec_liq(j)
186      sens_prec_sol_o(i,1) = sens_prec_sol(j)
187      lat_prec_liq_o(i,1) = lat_prec_liq(j)
188      lat_prec_sol_o(i,1) = lat_prec_sol(j)
189    enddo
190
191
192! - Flux calculation at first modele level for U and V
193    CALL calcul_flux_wind(knon, dtime, &
194         u0, v0, u1, v1, gustiness, cdragm, &
195         AcoefU, AcoefV, BcoefU, BcoefV, &
196         p1lay, temp_air, &
197         flux_u1, flux_v1) 
198 
199
200#ifdef ISO     
201        call calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
202     &    ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
203     &    evap, Roce,xtevap,h1 &
204#ifdef ISOTRAC
205     &    ,knindex &
206#endif
207     &   )
208#endif         
209
210#ifdef ISO
211#ifdef ISOVERIF
212!          write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
213          if (iso_eau.gt.0) then
214           do i=1,knon               
215              call iso_verif_egalite_choix(snow(i), &
216     &          xtsnow(iso_eau,i),'ocean_forced_mod 180', &
217     &          errmax,errmaxrel)
218           enddo ! do j=1,knon
219          endif !if (iso_eau.gt.0) then
220#endif
221#endif   
222
223  END SUBROUTINE ocean_forced_noice
224!
225!***************************************************************************************
226!
227  SUBROUTINE ocean_forced_ice( &
228       itime, dtime, jour, knon, knindex, &
229       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
230       AcoefH, AcoefQ, BcoefH, BcoefQ, &
231       AcoefU, AcoefV, BcoefU, BcoefV, &
232       ps, u1, v1, gustiness, &
233       radsol, snow, qsol, agesno, tsoil, &
234       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
235       tsurf_new, dflux_s, dflux_l &
236#ifdef ISO
237            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
238            xtsnow, xtsol,xtevap,Rland_ice & 
239#endif           
240            )
241!
242! This subroutine treats the ocean where there is ice.
243! The routine reads data from climatologie file and does flux calculations at the
244! surface.
245!
246    USE dimphy
247    USE calcul_fluxs_mod
248    USE surface_data,     ONLY : calice, calsno
249    USE limit_read_mod
250    USE fonte_neige_mod,  ONLY : fonte_neige
251    USE indice_sol_mod
252    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
253#ifdef ISO
254  USE infotrac_phy, ONLY: niso,ntraciso
255    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, &
256&       calcul_iso_surf_sic_vectall
257#ifdef ISOVERIF
258    USE isotopes_mod, ONLY: iso_eau,ridicule
259    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
260    USE isotopes_verif_mod
261#endif
262#endif
263
264!    INCLUDE "indicesol.h"
265    INCLUDE "dimsoil.h"
266    INCLUDE "YOMCST.h"
267    INCLUDE "clesphys.h"
268
269! Input arguments
270!****************************************************************************************
271    INTEGER, INTENT(IN)                  :: itime, jour, knon
272    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
273    REAL, INTENT(IN)                     :: dtime
274    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
275    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
276    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
277    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
278    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
279    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
280    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
281    REAL, DIMENSION(klon), INTENT(IN)    :: ps
282    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
283#ifdef ISO
284    REAL, DIMENSION(ntraciso,klon), INTENT(IN)    :: xtprecip_rain, xtprecip_snow
285    REAL, DIMENSION(ntraciso,klon), INTENT(IN)    :: xtspechum
286    REAL, DIMENSION(niso,klon), INTENT(IN)    :: Roce
287    REAL, DIMENSION(niso,klon), INTENT(IN)        :: Rland_ice
288#endif
289
290! In/Output arguments
291!****************************************************************************************
292    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
293    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
294    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
295    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
296#ifdef ISO     
297    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow
298    REAL, DIMENSION(niso,klon), INTENT(IN)    :: xtsol
299#endif
300
301! Output arguments
302!****************************************************************************************
303    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
304    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
305    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
306    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
307    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
308    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
309    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
310#ifdef ISO     
311    REAL, DIMENSION(ntraciso,klon), INTENT(OUT)    :: xtevap
312#endif     
313
314! Local variables
315!****************************************************************************************
316    LOGICAL                     :: check=.FALSE.
317    INTEGER                     :: i, j
318    REAL                        :: zfra
319    REAL, PARAMETER             :: t_grnd=271.35
320    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
321    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
322    REAL, DIMENSION(klon)       :: soilcap, soilflux
323    REAL, DIMENSION(klon)       :: u0, v0
324    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
325    REAL, DIMENSION(klon)       :: sens_prec_liq, sens_prec_sol   
326    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol   
327
328#ifdef ISO       
329      real, parameter :: t_coup = 273.15
330      real, dimension(klon) :: fq_fonte_diag
331      real, dimension(klon) :: fqfonte_diag
332      real, dimension(klon) ::  snow_evap_diag
333      real, dimension(klon) ::  fqcalving_diag
334      real, DIMENSION(klon) :: run_off_lic_diag
335      real :: coeff_rel_diag
336      real  :: max_eau_sol_diag 
337      real, dimension(klon)  ::  runoff_diag   
338      integer ixt
339      REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
340      REAL, DIMENSION(klon) :: snow_prec,qsol_prec 
341#endif     
342
343!****************************************************************************************
344! Start calculation
345!****************************************************************************************
346    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon
347
348!****************************************************************************************
349! 1)
350! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
351!                    dflux_s, dflux_l and qsurf
352!****************************************************************************************
353
354    tsurf_tmp(:) = tsurf_in(:)
355
356! calculate the parameters cal, beta, capsol and dif_grnd
357    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
358
359   
360    IF (soil_model) THEN
361! update tsoil and calculate soilcap and soilflux
362       CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, tsoil,soilcap, soilflux)
363       cal(1:knon) = RCPD / soilcap(1:knon)
364       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
365       dif_grnd = 1.0 / tau_gl
366    ELSE
367       dif_grnd = 1.0 / tau_gl
368       cal = RCPD * calice
369       WHERE (snow > 0.0) cal = RCPD * calsno
370    ENDIF
371
372    beta = 1.0
373    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
374
375! Suppose zero surface speed
376    u0(:)=0.0
377    v0(:)=0.0
378    u1_lay(:) = u1(:) - u0(:)
379    v1_lay(:) = v1(:) - v0(:)
380    CALL calcul_fluxs(knon, is_sic, dtime, &
381         tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, &
382         precip_rain, precip_snow, snow, qsurf,  &
383         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
384         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
385         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
386         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
387    do j = 1, knon
388      i = knindex(j)
389      sens_prec_liq_o(i,2) = sens_prec_liq(j)
390      sens_prec_sol_o(i,2) = sens_prec_sol(j)
391      lat_prec_liq_o(i,2) = lat_prec_liq(j)
392      lat_prec_sol_o(i,2) = lat_prec_sol(j)
393    enddo
394
395! - Flux calculation at first modele level for U and V
396    CALL calcul_flux_wind(knon, dtime, &
397         u0, v0, u1, v1, gustiness, cdragm, &
398         AcoefU, AcoefV, BcoefU, BcoefV, &
399         p1lay, temp_air, &
400         flux_u1, flux_v1) 
401
402!****************************************************************************************
403! 2)
404! Calculations due to snow and runoff
405!
406!****************************************************************************************
407
408#ifdef ISO
409   ! verif
410#ifdef ISOVERIF
411     do i=1,knon
412       if (iso_eau.gt.0) then
413           if (snow(i).gt.ridicule) then
414             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
415     &           'interfsurf 964',errmax,errmaxrel)
416            endif !if ((snow(i).gt.ridicule)) then
417        endif !if (iso_eau.gt.0) then     
418      enddo !do i=1,knon 
419#endif
420   ! end verif
421
422    do i=1,knon
423      snow_prec(i)=snow(i)
424      do ixt=1,niso
425      xtsnow_prec(ixt,i)=xtsnow(ixt,i)
426      enddo !do ixt=1,niso
427      ! initialisation:
428      fq_fonte_diag(i)=0.0
429      fqfonte_diag(i)=0.0
430      snow_evap_diag(i)=0.0
431    enddo !do i=1,knon
432#endif
433
434    CALL fonte_neige( knon, is_sic, knindex, dtime, &
435         tsurf_tmp, precip_rain, precip_snow, &
436         snow, qsol, tsurf_new, evap &
437#ifdef ISO   
438     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
439     & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
440#endif
441     &   )
442
443
444#ifdef ISO
445! isotopes: tout est externalisé
446#ifdef ISOVERIF
447        write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
448        write(*,*) 'klon,knon=',klon,knon
449#endif
450         call calcul_iso_surf_sic_vectall(klon,knon, &
451             &   evap,snow_evap_diag,Tsurf_new,Roce,snow, &
452             &   fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
453             &   precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
454             &   xtspechum,spechum,ps, &
455             &   xtevap,xtsnow,fqcalving_diag, &
456             &   knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
457     &  )   
458#ifdef ISOVERIF
459        write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
460          if (iso_eau.gt.0) then
461           do i=1,knon 
462              call iso_verif_egalite_choix(snow(i), &
463     &          xtsnow(iso_eau,i),'ocean_forced_mod 396', &
464     &          errmax,errmaxrel)
465           enddo ! do j=1,knon
466          endif !if (iso_eau.gt.0) then
467#endif
468!#ifdef ISOVERIF
469#endif   
470!#ifdef ISO
471   
472! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
473!
474    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 
475
476    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
477
478    alb1_new(:) = 0.0
479    DO i=1, knon
480       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
481       alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
482    ENDDO
483
484    alb2_new(:) = alb1_new(:)
485
486  END SUBROUTINE ocean_forced_ice
487
488!************************************************************************
489! 1D case
490!************************************************************************
491  SUBROUTINE read_tsurf1d(knon,sst_out)
492
493! This subroutine specifies the surface temperature to be used in 1D simulations
494
495      USE dimphy, ONLY : klon
496
497      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
498      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
499
500       INTEGER :: i
501! COMMON defined in lmdz1d.F:
502       real ts_cur
503       common /sst_forcing/ts_cur
504
505       DO i = 1, knon
506        sst_out(i) = ts_cur
507       ENDDO
508
509      END SUBROUTINE read_tsurf1d
510
511!
512!************************************************************************
513!
514END MODULE ocean_forced_mod
515
516
517
518
519
520
Note: See TracBrowser for help on using the repository browser.