source: LMDZ6/trunk/libf/phylmdiso/surf_ocean_mod.F90 @ 3927

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

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 13.0 KB
Line 
1!
2! $Id: surf_ocean_mod.F90 3395 2018-09-25 15:22:13Z lguez $
3!
4MODULE surf_ocean_mod
5
6  IMPLICIT NONE
7
8CONTAINS
9  !
10  !******************************************************************************
11  !
12  SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
13       windsp, rmu0, fder, tsurf_in, &
14       itime, dtime, jour, knon, knindex, &
15       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
16       AcoefH, AcoefQ, BcoefH, BcoefQ, &
17       AcoefU, AcoefV, BcoefU, BcoefV, &
18       ps, u1, v1, gustiness, rugoro, pctsrf, &
19       snow, qsurf, agesno, &
20       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &
21       tsurf_new, dflux_s, dflux_l, lmt_bils, &
22       flux_u1, flux_v1 &
23#ifdef ISO
24        &       ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
25        &       xtsnow,xtevap,h1 &
26#endif               
27        &       )
28
29    use albedo, only: alboc, alboc_cd
30    USE dimphy, ONLY: klon, zmasq
31    USE surface_data, ONLY     : type_ocean
32    USE ocean_forced_mod, ONLY : ocean_forced_noice
33    USE ocean_slab_mod, ONLY   : ocean_slab_noice
34    USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
35    USE indice_sol_mod, ONLY : nbsrf, is_oce
36#ifdef ISO
37  USE infotrac_phy, ONLY : ntraciso,niso
38#ifdef ISOVERIF
39    USE isotopes_mod, ONLY: iso_eau,ridicule
40    USE isotopes_verif_mod
41#endif
42#endif
43    USE limit_read_mod
44    !
45    ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force,
46    ! slab or couple). The calculations of albedo and rugosity for the ocean surface are
47    ! done in here because they are identical for the different modes of ocean.
48
49
50    INCLUDE "YOMCST.h"
51
52    include "clesphys.h"
53    ! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0)
54
55    ! Input variables
56    !******************************************************************************
57    INTEGER, INTENT(IN)                      :: itime, jour, knon
58    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
59    REAL, INTENT(IN)                         :: dtime
60    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
61    REAL, DIMENSION(klon), INTENT(IN)        :: swnet  ! net shortwave radiation at surface 
62    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface 
63    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
64    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
65    REAL, DIMENSION(klon), INTENT(IN)        :: rmu0 
66    REAL, DIMENSION(klon), INTENT(IN)        :: fder
67    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
68    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
69    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
70    REAL, DIMENSION(klon), INTENT(IN)        :: cdragm
71    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
72    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
73    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
74    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
75    REAL, DIMENSION(klon), INTENT(IN)        :: ps
76    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
77    REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
78    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
79#ifdef ISO
80    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: xtprecip_rain, xtprecip_snow
81    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: xtspechum
82#endif
83
84    ! In/Output variables
85    !******************************************************************************
86    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
87    REAL, DIMENSION(klon), INTENT(INOUT)     :: qsurf
88    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
89#ifdef ISO
90    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsnow
91    REAL, DIMENSION(niso,klon), INTENT(INOUT)        :: Roce 
92#endif
93    REAL, DIMENSION(klon), INTENT(inOUT):: z0h
94
95    ! Output variables
96    !******************************************************************************
97    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m
98    !albedo SB >>>
99    !    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
100    !    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
101    REAL, DIMENSION(6), INTENT(IN)          :: SFRWL
102    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
103    !albedo SB <<<     
104    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
105    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
106    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
107    REAL, DIMENSION(klon), INTENT(OUT)       :: lmt_bils
108    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
109#ifdef ISO
110    REAL, DIMENSION(ntraciso,klon), INTENT(out)        :: xtevap ! isotopes in surface evaporation flux
111    REAL, DIMENSION(klon), INTENT(out)        :: h1 ! just a diagnostic, not useful for the simulation   
112#endif
113
114    ! Local variables
115    !******************************************************************************
116    INTEGER               :: i, k
117    REAL                  :: tmp
118    REAL, PARAMETER       :: cepdu2=(0.1)**2
119    REAL, DIMENSION(klon) :: alb_eau, z0_lim
120    REAL, DIMENSION(klon) :: radsol
121    REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
122    CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
123
124    ! End definition
125    !******************************************************************************
126
127#ifdef ISO
128#ifdef ISOVERIF
129        do i=1,knon
130          if (iso_eau.gt.0) then         
131                 call iso_verif_egalite_choix(xtspechum(iso_eau,i), &
132     &                  spechum(i),'surf_ocean_mod 117', &
133     &                  errmax,errmaxrel)         
134                 call iso_verif_egalite_choix(xtsnow(iso_eau,i), &
135     &                  snow(i),'surf_ocean_mod 127', &
136     &                  errmax,errmaxrel)
137           endif !if (iso_eau.gt.0) then
138        enddo !do i=1,klon
139#endif     
140#endif
141
142    !******************************************************************************
143    ! Calculate total net radiance at surface
144    !
145    !******************************************************************************
146    radsol(1:klon) = 0.0 ! initialisation a priori inutile
147    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
148
149    !******************************************************************************
150    ! Cdragq computed from cdrag
151    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
152    ! it can be computed inside surf_ocean
153    ! More complicated appraches may require the propagation through
154    ! pbl_surface of an independant cdragq variable.
155    !******************************************************************************
156
157    IF ( f_z0qh_oce .ne. 1.) THEN
158       ! Si on suit les formulations par exemple de Tessel, on
159       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
160       cdragq(1:knon)=cdragh(1:knon)*                                      &
161            log(z1lay(1:knon)/z0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*z0h(1:knon)))
162    ELSE
163       cdragq(1:knon)=cdragh(1:knon)
164    ENDIF
165
166
167    !******************************************************************************
168    ! Switch according to type of ocean (couple, slab or forced)
169    !******************************************************************************
170    SELECT CASE(type_ocean)
171    CASE('couple')
172       CALL ocean_cpl_noice( &
173            swnet, lwnet, alb1, &
174            windsp, fder, &
175            itime, dtime, knon, knindex, &
176            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,&
177            AcoefH, AcoefQ, BcoefH, BcoefQ, &
178            AcoefU, AcoefV, BcoefU, BcoefV, &
179            ps, u1, v1, gustiness, &
180            radsol, snow, agesno, &
181            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
182            tsurf_new, dflux_s, dflux_l)
183
184    CASE('slab')
185       CALL ocean_slab_noice( &
186            itime, dtime, jour, knon, knindex, &
187            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
188            AcoefH, AcoefQ, BcoefH, BcoefQ, &
189            AcoefU, AcoefV, BcoefU, BcoefV, &
190            ps, u1, v1, gustiness, tsurf_in, &
191            radsol, snow, &
192            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
193            tsurf_new, dflux_s, dflux_l, lmt_bils)
194
195    CASE('force')
196       CALL ocean_forced_noice( &
197            itime, dtime, jour, knon, knindex, &
198            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
199            temp_air, spechum, &
200            AcoefH, AcoefQ, BcoefH, BcoefQ, &
201            AcoefU, AcoefV, BcoefU, BcoefV, &
202            ps, u1, v1, gustiness, &
203            radsol, snow, agesno, &
204            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
205            tsurf_new, dflux_s, dflux_l &
206#ifdef ISO
207            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
208            xtsnow,xtevap,h1 & 
209#endif           
210            )
211    END SELECT
212
213    !******************************************************************************
214    ! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
215    !******************************************************************************
216    IF (type_ocean.NE.'slab') THEN
217       lmt_bils(1:klon)=0.
218       DO i=1,knon
219          lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
220               *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i)))
221       END DO
222    END IF
223
224    !******************************************************************************
225    ! Calculate ocean surface albedo
226    !******************************************************************************
227    !albedo SB >>>
228    IF (iflag_albedo==0) THEN
229       !--old parametrizations of ocean surface albedo
230       !
231       IF (iflag_cycle_diurne.GE.1) THEN
232          !
233          CALL alboc_cd(rmu0,alb_eau)
234          !
235          !--ad-hoc correction for model radiative balance tuning
236          !--now outside alboc_cd routine
237          alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
238          alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.0),1.0)
239          !
240       ELSE
241          !
242          CALL alboc(REAL(jour),rlat,alb_eau)
243          !--ad-hoc correction for model radiative balance tuning
244          !--now outside alboc routine
245          alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
246          alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.04),0.60)
247          !
248       ENDIF
249       !
250       DO i =1, knon
251          DO  k=1,nsw
252             alb_dir_new(i,k) = alb_eau(knindex(i))
253          ENDDO
254       ENDDO
255       !IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions
256       !albedo for diffuse radiation is taken the same as for direct radiation
257       alb_dif_new(1:knon,:)=alb_dir_new(1:knon,:)
258       !IM 09122015 end
259       !
260    ELSE IF (iflag_albedo==1) THEN
261       !--new parametrization of ocean surface albedo by Sunghye Baek
262       !--albedo for direct and diffuse radiation are different
263       !
264       CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
265       !
266       !--ad-hoc correction for model radiative balance tuning
267       alb_dir_new(1:knon,:) = fmagic*alb_dir_new(1:knon,:) + pmagic
268       alb_dif_new(1:knon,:) = fmagic*alb_dif_new(1:knon,:) + pmagic
269       alb_dir_new(1:knon,:)=MIN(MAX(alb_dir_new(1:knon,:),0.0),1.0)
270       alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
271       !
272    ELSE IF (iflag_albedo==2) THEN
273       ! F. Codron albedo read from limit.nc
274       CALL limit_read_rug_alb(itime, dtime, jour,&
275            knon, knindex, z0_lim, alb_eau)
276       DO i =1, knon
277          DO  k=1,nsw
278             alb_dir_new(i,k) = alb_eau(i)
279          ENDDO
280       ENDDO
281       alb_dif_new=alb_dir_new
282    ENDIF
283    !albedo SB <<<
284
285    !******************************************************************************
286    ! Calculate the rugosity
287    !******************************************************************************
288    IF (iflag_z0_oce==0) THEN
289       DO i = 1, knon
290          tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
291          z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
292               +  0.11*14e-6 / SQRT(cdragm(i) * tmp)
293          z0m(i) = MAX(1.5e-05,z0m(i))
294       ENDDO
295       z0h(1:knon)=z0m(1:knon) ! En attendant mieux
296
297    ELSE IF (iflag_z0_oce==1) THEN
298       DO i = 1, knon
299          tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
300          z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
301               + 0.11*14e-6 / SQRT(cdragm(i) * tmp)
302          z0m(i) = MAX(1.5e-05,z0m(i))
303          z0h(i)=0.4*14e-6 / SQRT(cdragm(i) * tmp)
304       ENDDO
305    ELSE IF (iflag_z0_oce==-1) THEN
306       DO i = 1, knon
307          z0m(i) = z0min
308          z0h(i) = z0min
309       ENDDO
310    ELSE
311       CALL abort_physic(modname,'version non prevue',1)
312    ENDIF
313    !
314    !******************************************************************************
315  END SUBROUTINE surf_ocean
316  !******************************************************************************
317  !
318END MODULE surf_ocean_mod
Note: See TracBrowser for help on using the repository browser.