source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyaqua_mod.F90

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

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 32.1 KB
Line 
1!
2! $Id: phyaqua_mod.F90 3579 2019-10-09 13:11:07Z fairhead $
3!
4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE geometry_mod, ONLY : latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29#ifdef ISO
30  USE fonte_neige_mod,  ONLY : fonte_neige_init_iso
31  USE pbl_surface_mod,  ONLY : pbl_surface_init_iso
32#endif
33    USE phys_state_var_mod
34    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
35                                day_ini,day_end
36    USE indice_sol_mod
37    USE nrtype, ONLY: pi
38!    USE ioipsl
39    USE mod_phys_lmdz_para, ONLY: is_master
40    USE mod_phys_lmdz_transfert_para, ONLY: bcast
41    USE mod_grid_phy_lmdz
42    USE ioipsl_getin_p_mod, ONLY : getin_p
43    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
44#ifdef ISO
45  USE infotrac_phy, ONLY: niso
46#endif
47
48    IMPLICIT NONE
49
50    include "YOMCST.h"
51    include "clesphys.h"
52    include "dimsoil.h"
53
54    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
55    ! IM ajout latfi, lonfi
56!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
57
58    INTEGER type_profil, type_aqua
59
60    ! Ajouts initialisation des surfaces
61    REAL :: run_off_lic_0(nlon)
62    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
63    REAL :: tsoil(nlon, nsoilmx, nbsrf)
64    REAL :: tslab(nlon), seaice(nlon)
65    REAL fder(nlon)
66
67#ifdef ISO
68    REAL :: xtrun_off_lic_0(niso,nlon)
69    REAL :: xtsolsrf(niso,nlon, nbsrf), xtsnsrf(niso,nlon, nbsrf)
70    REAL :: Rland_ice(niso,nlon)
71#endif
72
73
74
75    ! Arguments :
76    ! -----------
77
78    ! integer radpas
79    INTEGER it, unit, i, k, itap
80
81    REAL rugos, albedo
82    REAL tsurf
83    REAL time, timestep, day, day0
84    REAL qsol_f
85    REAL rugsrel(nlon)
86    LOGICAL alb_ocean
87
88    CHARACTER *80 ans, file_forctl, file_fordat, file_start
89    CHARACTER *100 file, var
90    CHARACTER *2 cnbl
91
92    REAL phy_nat(nlon, year_len)
93    REAL phy_alb(nlon, year_len)
94    REAL phy_sst(nlon, year_len)
95    REAL phy_bil(nlon, year_len)
96    REAL phy_rug(nlon, year_len)
97    REAL phy_ice(nlon, year_len)
98    REAL phy_fter(nlon, year_len)
99    REAL phy_foce(nlon, year_len)
100    REAL phy_fsic(nlon, year_len)
101    REAL phy_flic(nlon, year_len)
102
103    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
104!$OMP THREADPRIVATE(read_climoz)
105
106    ! -------------------------------------------------------------------------
107    ! declaration pour l'appel a phyredem
108    ! -------------------------------------------------------------------------
109
110    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
111    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
112    ! real pbl_tke(nlon,llm,nbsrf)
113    ! real rain_fall(nlon),snow_fall(nlon)
114    ! real solsw(nlon), sollw(nlon),radsol(nlon)
115    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
116    ! real ratqs(nlon,llm)
117    ! real clwcon(nlon,llm)
118
119    INTEGER longcles
120    PARAMETER (longcles=20)
121    REAL clesphy0(longcles)
122
123
124    ! -----------------------------------------------------------------------
125    ! dynamial tendencies :
126    ! ---------------------
127
128    INTEGER l, ierr, aslun
129
130    REAL paire
131
132    ! Local
133    CHARACTER (LEN=20) :: modname='phyaqua'
134    CHARACTER (LEN=80) :: abort_message
135
136
137    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
138    ! INITIALISATIONS
139    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140
141    ! -----------------------------------------------------------------------
142    ! Initialisations  des constantes
143    ! -------------------------------
144
145    !IF (calend .EQ. "earth_360d") Then
146      year_len_phy = year_len
147    !END IF
148   
149    if (year_len.ne.360) then
150      write (*,*) year_len
151      write (*,*) 'iniaqua: 360 day calendar is required !'
152      stop
153    endif
154
155    type_aqua = iflag_phys/100
156    type_profil = iflag_phys - type_aqua*100
157    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
158
159    IF (klon/=nlon) THEN
160      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
161      abort_message= 'probleme de dimensions dans iniaqua'
162      CALL abort_physic(modname,abort_message,1)
163    END IF
164    CALL phys_state_var_init(read_climoz)
165
166
167    read_climoz = 0
168    day0 = 217.
169    day = day0
170    it = 0
171    time = 0.
172
173    ! -----------------------------------------------------------------------
174    ! initialisations de la physique
175    ! -----------------------------------------------------------------------
176
177    day_ini = day_ref
178    day_end = day_ini + ndays
179
180    nbapp_rad = 24
181    CALL getin_p('nbapp_rad', nbapp_rad)
182
183    ! ---------------------------------------------------------------------
184    ! Creation des conditions aux limites:
185    ! ------------------------------------
186    ! Initialisations des constantes
187    ! Ajouter les manquants dans planete.def... (albedo etc)
188    co2_ppm = 348.
189    CALL getin_p('co2_ppm', co2_ppm)
190
191    solaire = 1365.
192    CALL getin_p('solaire', solaire)
193 
194    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
195    ! type_aqua
196    alb_ocean = .TRUE.
197    CALL getin_p('alb_ocean', alb_ocean)
198
199    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
200    WRITE (*, *) 'iniaqua: solaire=', solaire
201    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
202
203    radsol = 0.
204    qsol_f = 10.
205
206    ! Conditions aux limites:
207    ! -----------------------
208
209    qsol(:) = qsol_f
210    rugsrel = 0.0 ! (rugsrel = rugoro)
211    rugoro = 0.0
212    u_ancien = 0.0
213    v_ancien = 0.0
214    agesno = 50.0
215    ! Relief plat
216    zmea = 0.
217    zstd = 0.
218    zsig = 0.
219    zgam = 0.
220    zthe = 0.
221    zpic = 0.
222    zval = 0.
223
224    ! Une seule surface
225    pctsrf = 0.
226    IF (type_aqua==1) THEN
227      rugos = 1.E-4
228      albedo = 0.19
229      pctsrf(:, is_oce) = 1.
230    ELSE IF (type_aqua==2) THEN
231      rugos = 0.03
232      albedo = 0.1
233      pctsrf(:, is_ter) = 1.
234    END IF
235
236    CALL getin_p('rugos', rugos)
237
238    WRITE (*, *) 'iniaqua: rugos=', rugos
239    zmasq(:) = pctsrf(:, is_ter)
240
241    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
242    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
243
244    ! Si alb_ocean on calcule un albedo oceanique moyen
245    ! if (alb_ocean) then
246    ! Voir pourquoi on avait ca.
247    ! CALL ini_alb_oce(phy_alb)
248    ! else
249    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
250    ! endif !alb_ocean
251
252    DO i = 1, year_len
253      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
254      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
255
256      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
257      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
258      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
259      phy_ice(:, i) = 0.0 ! fraction de glace (?)
260      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
261      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
262      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
263      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
264    END DO
265    ! IM calcul profil sst
266    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
267
268    IF (grid_type==unstructured) THEN
269      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
270                             phy_fter, phy_foce, phy_flic, phy_fsic)
271    ELSE
272     
273       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
274                     phy_fter, phy_foce, phy_flic, phy_fsic)
275    ENDIF
276
277    ! ---------------------------------------------------------------------
278    ! Ecriture de l'etat initial:
279    ! ---------------------------
280
281
282    ! Ecriture etat initial physique
283
284    timestep = pdtphys
285    radpas = nint(rday/timestep/float(nbapp_rad))
286
287    DO i = 1, longcles
288      clesphy0(i) = 0.
289    END DO
290    clesphy0(1) = float(iflag_con)
291    clesphy0(2) = float(nbapp_rad)
292    ! IF( cycle_diurne  ) clesphy0(3) =  1.
293    clesphy0(3) = 1. ! cycle_diurne
294    clesphy0(4) = 1. ! soil_model
295    clesphy0(5) = 1. ! new_oliq
296    clesphy0(6) = 0. ! ok_orodr
297    clesphy0(7) = 0. ! ok_orolf
298    clesphy0(8) = 0. ! ok_limitvrai
299
300
301    ! =======================================================================
302    ! Profils initiaux
303    ! =======================================================================
304
305    ! On initialise les temperatures de surfaces comme les sst
306    DO i = 1, nlon
307      ftsol(i, :) = phy_sst(i, 1)
308      tsoil(i, :, :) = phy_sst(i, 1)
309      tslab(i) = phy_sst(i, 1)
310    END DO
311
312    falbe(:, :) = albedo
313    falblw(:, :) = albedo
314    rain_fall(:) = 0.
315    snow_fall(:) = 0.
316    solsw(:) = 0.
317    sollw(:) = 0.
318    radsol(:) = 0.
319
320    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
321    ! intialisation bidon mais pas grave
322    t_ancien(:, :) = 0.
323    q_ancien(:, :) = 0.
324    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
325    rnebcon = 0.
326    ratqs = 0.
327    clwcon = 0.
328    pbl_tke = 1.E-8
329
330    ! variables supplementaires pour appel a plb_surface_init
331    fder(:) = 0.
332    seaice(:) = 0.
333    run_off_lic_0 = 0.
334    fevap = 0.
335#ifdef ISO
336    xtrun_off_lic_0 = 0.
337#endif
338
339
340    ! Initialisations necessaires avant phyredem
341    type_ocean = 'force'
342    CALL fonte_neige_init(run_off_lic_0)
343#ifdef ISO
344   CALL fonte_neige_init_iso(xtrun_off_lic_0)
345#endif
346    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
347    snsrf(:, :) = 0. ! couverture de neige des sous surface
348    z0m(:, :) = rugos ! couverture de neige des sous surface
349    z0h=z0m
350
351
352    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
353#ifdef ISO
354  CALL pbl_surface_init_iso(xtsnsrf,Rland_ice)
355#endif
356
357    PRINT *, 'iniaqua: before phyredem'
358
359    pbl_tke(:,:,:) = 1.e-8
360    falb1 = albedo
361    falb2 = albedo
362    zmax0 = 0.
363    f0 = 0.
364    sig1 = 0.
365    w01 = 0.
366    wake_deltat = 0.
367    wake_deltaq = 0.
368    wake_s = 0.
369    wake_dens = 0.
370    wake_cstar = 0.
371    wake_pe = 0.
372    wake_fip = 0.
373    fm_therm = 0.
374    entr_therm = 0.
375    detr_therm = 0.
376    ale_bl = 0.
377    ale_bl_trig =0.
378    alp_bl =0.
379    treedrg(:,:,:)=0.
380
381    u10m = 0.
382    v10m = 0.
383
384    ql_ancien   = 0.
385    qs_ancien   = 0.
386    u_ancien    = 0.
387    v_ancien    = 0.
388    prw_ancien  = 0.
389    prlw_ancien = 0.
390    prsw_ancien = 0. 
391
392    ale_wake    = 0.
393    ale_bl_stat = 0. 
394
395
396!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
397!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
398!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
399    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
400    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
401
402!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
403!ym probably the uninitialized value was 0 for standard (regular grid) case
404    falb_dif(:,:,:)=0
405
406
407    CALL phyredem('startphy.nc')
408
409    PRINT *, 'iniaqua: after phyredem'
410    CALL phys_state_var_end
411
412    RETURN
413  END SUBROUTINE iniaqua
414
415
416  ! ====================================================================
417  ! ====================================================================
418  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
419    USE dimphy
420    IMPLICIT NONE
421    ! ====================================================================
422    ! =============================================================
423    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
424    ! Auteur : A. Campoy et F. Hourdin
425    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
426    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
427    ! connaissant la declinaison, la latitude et la longitude.
428
429    ! Dans cette version particuliere, on calcule le rayonnement
430    ! moyen sur l'année à chaque latitude.
431    ! angle zenithal calculé pour obtenir un
432    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
433    ! en moyenne annuelle.
434    ! Spécifique de la terre. Utilisé pour les aqua planetes.
435
436    ! Rque   : Different de la routine angle en ce sens que zenang
437    ! fournit des moyennes de pmu0 et non des valeurs
438    ! instantanees, du coup frac prend toutes les valeurs
439    ! entre 0 et 1.
440    ! Date   : premiere version le 13 decembre 1994
441    ! revu pour  GCM  le 30 septembre 1996
442    ! ===============================================================
443    ! longi----INPUT : la longitude vraie de la terre dans son plan
444    ! solaire a partir de l'equinoxe de printemps (degre)
445    ! gmtime---INPUT : temps universel en fraction de jour
446    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
447    ! lat------INPUT : latitude en degres
448    ! long-----INPUT : longitude en degres
449    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
450    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
451    ! ================================================================
452    include "YOMCST.h"
453    ! ================================================================
454    LOGICAL cycle_diurne
455    REAL gmtime
456    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
457    ! ================================================================
458    INTEGER i
459    REAL gmtime1, gmtime2
460    REAL pi_local
461
462
463    REAL rmu0m(klon), rmu0a(klon)
464
465
466    pi_local = 4.0*atan(1.0)
467
468    ! ================================================================
469    ! Calcul de l'angle zenithal moyen sur la journee
470    ! ================================================================
471
472    DO i = 1, klon
473      fract(i) = 1.
474      ! Calcule du flux moyen
475      IF (abs(rlat(i))<=28.75) THEN
476        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
477      ELSE IF (abs(rlat(i))<=43.75) THEN
478        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
479      ELSE IF (abs(rlat(i))<=71.25) THEN
480        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
481      ELSE
482        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
483      END IF
484    END DO
485
486    ! ================================================================
487    ! Avec ou sans cycle diurne
488    ! ================================================================
489
490    IF (cycle_diurne) THEN
491
492      ! On redecompose flux  au sommet suivant un cycle diurne idealise
493      ! identique a toutes les latitudes.
494
495      DO i = 1, klon
496        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
497        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
498          rmu0a(i)/sqrt(2.)
499      END DO
500
501      DO i = 1, klon
502        IF (rmu0(i)<=0.) THEN
503          rmu0(i) = 0.
504          fract(i) = 0.
505        ELSE
506          fract(i) = 1.
507        END IF
508      END DO
509
510      ! Affichage de l'angel zenitale
511      ! print*,'************************************'
512      ! print*,'************************************'
513      ! print*,'************************************'
514      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
515      ! print*,'rmu0m=',rmu0m(i)
516      ! print*,'rmu0a=',rmu0a(i)
517      ! print*,'rmu0=',rmu0(i)
518
519    ELSE
520
521      DO i = 1, klon
522        fract(i) = 0.5
523        rmu0(i) = rmu0m(i)*2.
524      END DO
525
526    END IF
527
528    RETURN
529  END SUBROUTINE zenang_an
530
531  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
532
533  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
534      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
535
536    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
537    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
538#ifdef CPP_XIOS
539    USE xios
540#endif
541    IMPLICIT NONE
542
543    include "netcdf.inc"
544
545    INTEGER, INTENT (IN) :: klon
546    REAL, INTENT (IN) :: phy_nat(klon, 360)
547    REAL, INTENT (IN) :: phy_alb(klon, 360)
548    REAL, INTENT (IN) :: phy_sst(klon, 360)
549    REAL, INTENT (IN) :: phy_bil(klon, 360)
550    REAL, INTENT (IN) :: phy_rug(klon, 360)
551    REAL, INTENT (IN) :: phy_ice(klon, 360)
552    REAL, INTENT (IN) :: phy_fter(klon, 360)
553    REAL, INTENT (IN) :: phy_foce(klon, 360)
554    REAL, INTENT (IN) :: phy_flic(klon, 360)
555    REAL, INTENT (IN) :: phy_fsic(klon, 360)
556
557    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
558      ! on the whole physics grid
559 
560#ifdef CPP_XIOS
561    PRINT *, 'writelim: Ecriture du fichier limit'
562
563    CALL gather_omp(phy_foce, phy_mpi)
564    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
565
566    CALL gather_omp(phy_fsic, phy_mpi)
567    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
568     
569    CALL gather_omp(phy_fter, phy_mpi)
570    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
571     
572    CALL gather_omp(phy_flic, phy_mpi)
573    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
574
575    CALL gather_omp(phy_sst, phy_mpi)
576    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
577
578    CALL gather_omp(phy_bil, phy_mpi)
579    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
580
581    CALL gather_omp(phy_alb, phy_mpi)
582    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
583
584    CALL gather_omp(phy_rug, phy_mpi)
585    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
586#endif
587  END SUBROUTINE writelim_unstruct
588
589
590
591  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
592      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
593
594    USE mod_phys_lmdz_para, ONLY: is_master
595    USE mod_grid_phy_lmdz, ONLY: klon_glo
596    USE mod_phys_lmdz_transfert_para, ONLY: gather
597    USE phys_cal_mod, ONLY: year_len
598    IMPLICIT NONE
599    include "netcdf.inc"
600
601    INTEGER, INTENT (IN) :: klon
602    REAL, INTENT (IN) :: phy_nat(klon, year_len)
603    REAL, INTENT (IN) :: phy_alb(klon, year_len)
604    REAL, INTENT (IN) :: phy_sst(klon, year_len)
605    REAL, INTENT (IN) :: phy_bil(klon, year_len)
606    REAL, INTENT (IN) :: phy_rug(klon, year_len)
607    REAL, INTENT (IN) :: phy_ice(klon, year_len)
608    REAL, INTENT (IN) :: phy_fter(klon, year_len)
609    REAL, INTENT (IN) :: phy_foce(klon, year_len)
610    REAL, INTENT (IN) :: phy_flic(klon, year_len)
611    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
612
613    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
614      ! on the whole physics grid
615    INTEGER :: k
616    INTEGER ierr
617    INTEGER dimfirst(3)
618    INTEGER dimlast(3)
619
620    INTEGER nid, ndim, ntim
621    INTEGER dims(2), debut(2), epais(2)
622    INTEGER id_tim
623    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
624    INTEGER id_fter, id_foce, id_fsic, id_flic
625
626    IF (is_master) THEN
627
628      PRINT *, 'writelim: Ecriture du fichier limit'
629
630      ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
631
632      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
633        'Fichier conditions aux limites')
634      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
635      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
636      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
637
638      dims(1) = ndim
639      dims(2) = ntim
640
641#ifdef NC_DOUBLE
642      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
643#else
644      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
645#endif
646      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
647
648#ifdef NC_DOUBLE
649      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
650#else
651      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
652#endif
653      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
654        'Nature du sol (0,1,2,3)')
655
656#ifdef NC_DOUBLE
657      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
658#else
659      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
660#endif
661      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
662        'Temperature superficielle de la mer')
663
664#ifdef NC_DOUBLE
665      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
666#else
667      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
668#endif
669      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
670        'Reference flux de chaleur au sol')
671
672#ifdef NC_DOUBLE
673      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
674#else
675      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
676#endif
677      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
678
679#ifdef NC_DOUBLE
680      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
681#else
682      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
683#endif
684      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
685
686#ifdef NC_DOUBLE
687      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
688#else
689      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
690#endif
691      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
692#ifdef NC_DOUBLE
693      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
694#else
695      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
696#endif
697      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
698#ifdef NC_DOUBLE
699      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
700#else
701      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
702#endif
703      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
704#ifdef NC_DOUBLE
705      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
706#else
707      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
708#endif
709      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
710
711      ierr = nf_enddef(nid)
712      IF (ierr/=nf_noerr) THEN
713        WRITE (*, *) 'writelim error: failed to end define mode'
714        WRITE (*, *) nf_strerror(ierr)
715      END IF
716
717
718      ! write the 'times'
719      DO k = 1, year_len
720#ifdef NC_DOUBLE
721        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
722#else
723        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
724#endif
725        IF (ierr/=nf_noerr) THEN
726          WRITE (*, *) 'writelim error with temps(k),k=', k
727          WRITE (*, *) nf_strerror(ierr)
728        END IF
729      END DO
730
731    END IF ! of if (is_master)
732
733    ! write the fields, after having collected them on master
734
735    CALL gather(phy_nat, phy_glo)
736    IF (is_master) THEN
737#ifdef NC_DOUBLE
738      ierr = nf_put_var_double(nid, id_nat, phy_glo)
739#else
740      ierr = nf_put_var_real(nid, id_nat, phy_glo)
741#endif
742      IF (ierr/=nf_noerr) THEN
743        WRITE (*, *) 'writelim error with phy_nat'
744        WRITE (*, *) nf_strerror(ierr)
745      END IF
746    END IF
747
748    CALL gather(phy_sst, phy_glo)
749    IF (is_master) THEN
750#ifdef NC_DOUBLE
751      ierr = nf_put_var_double(nid, id_sst, phy_glo)
752#else
753      ierr = nf_put_var_real(nid, id_sst, phy_glo)
754#endif
755      IF (ierr/=nf_noerr) THEN
756        WRITE (*, *) 'writelim error with phy_sst'
757        WRITE (*, *) nf_strerror(ierr)
758      END IF
759    END IF
760
761    CALL gather(phy_bil, phy_glo)
762    IF (is_master) THEN
763#ifdef NC_DOUBLE
764      ierr = nf_put_var_double(nid, id_bils, phy_glo)
765#else
766      ierr = nf_put_var_real(nid, id_bils, phy_glo)
767#endif
768      IF (ierr/=nf_noerr) THEN
769        WRITE (*, *) 'writelim error with phy_bil'
770        WRITE (*, *) nf_strerror(ierr)
771      END IF
772    END IF
773
774    CALL gather(phy_alb, phy_glo)
775    IF (is_master) THEN
776#ifdef NC_DOUBLE
777      ierr = nf_put_var_double(nid, id_alb, phy_glo)
778#else
779      ierr = nf_put_var_real(nid, id_alb, phy_glo)
780#endif
781      IF (ierr/=nf_noerr) THEN
782        WRITE (*, *) 'writelim error with phy_alb'
783        WRITE (*, *) nf_strerror(ierr)
784      END IF
785    END IF
786
787    CALL gather(phy_rug, phy_glo)
788    IF (is_master) THEN
789#ifdef NC_DOUBLE
790      ierr = nf_put_var_double(nid, id_rug, phy_glo)
791#else
792      ierr = nf_put_var_real(nid, id_rug, phy_glo)
793#endif
794      IF (ierr/=nf_noerr) THEN
795        WRITE (*, *) 'writelim error with phy_rug'
796        WRITE (*, *) nf_strerror(ierr)
797      END IF
798    END IF
799
800    CALL gather(phy_fter, phy_glo)
801    IF (is_master) THEN
802#ifdef NC_DOUBLE
803      ierr = nf_put_var_double(nid, id_fter, phy_glo)
804#else
805      ierr = nf_put_var_real(nid, id_fter, phy_glo)
806#endif
807      IF (ierr/=nf_noerr) THEN
808        WRITE (*, *) 'writelim error with phy_fter'
809        WRITE (*, *) nf_strerror(ierr)
810      END IF
811    END IF
812
813    CALL gather(phy_foce, phy_glo)
814    IF (is_master) THEN
815#ifdef NC_DOUBLE
816      ierr = nf_put_var_double(nid, id_foce, phy_glo)
817#else
818      ierr = nf_put_var_real(nid, id_foce, phy_glo)
819#endif
820      IF (ierr/=nf_noerr) THEN
821        WRITE (*, *) 'writelim error with phy_foce'
822        WRITE (*, *) nf_strerror(ierr)
823      END IF
824    END IF
825
826    CALL gather(phy_fsic, phy_glo)
827    IF (is_master) THEN
828#ifdef NC_DOUBLE
829      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
830#else
831      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
832#endif
833      IF (ierr/=nf_noerr) THEN
834        WRITE (*, *) 'writelim error with phy_fsic'
835        WRITE (*, *) nf_strerror(ierr)
836      END IF
837    END IF
838
839    CALL gather(phy_flic, phy_glo)
840    IF (is_master) THEN
841#ifdef NC_DOUBLE
842      ierr = nf_put_var_double(nid, id_flic, phy_glo)
843#else
844      ierr = nf_put_var_real(nid, id_flic, phy_glo)
845#endif
846      IF (ierr/=nf_noerr) THEN
847        WRITE (*, *) 'writelim error with phy_flic'
848        WRITE (*, *) nf_strerror(ierr)
849      END IF
850    END IF
851
852    ! close file:
853    IF (is_master) THEN
854      ierr = nf_close(nid)
855    END IF
856
857  END SUBROUTINE writelim
858
859  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
860
861  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
862    USE dimphy
863    USE phys_cal_mod , ONLY: year_len
864    IMPLICIT NONE
865
866    INTEGER nlon, type_profil, i, k, j
867    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
868    INTEGER imn, imx, amn, amx, kmn, kmx
869    INTEGER p, pplus, nlat_max
870    PARAMETER (nlat_max=72)
871    REAL x_anom_sst(nlat_max)
872    CHARACTER (LEN=20) :: modname='profil_sst'
873    CHARACTER (LEN=80) :: abort_message
874
875    IF (klon/=nlon) THEN
876       abort_message='probleme de dimensions dans profil_sst'
877       CALL abort_physic(modname,abort_message,1)
878    ENDIF
879    WRITE (*, *) ' profil_sst: type_profil=', type_profil
880    DO i = 1, year_len
881      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
882
883      ! Rajout fbrlmd
884
885      IF (type_profil==1) THEN
886        ! Méthode 1 "Control" faible plateau à l'Equateur
887        DO j = 1, klon
888          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
889          ! PI/3=1.047197551
890          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
891            phy_sst(j, i) = 273.
892          END IF
893        END DO
894      END IF
895      IF (type_profil==2) THEN
896        ! Méthode 2 "Flat" fort plateau à l'Equateur
897        DO j = 1, klon
898          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
899          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
900            phy_sst(j, i) = 273.
901          END IF
902        END DO
903      END IF
904
905
906      IF (type_profil==3) THEN
907        ! Méthode 3 "Qobs" plateau réel à l'Equateur
908        DO j = 1, klon
909          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
910            rlatd(j))**4)
911          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
912            phy_sst(j, i) = 273.
913          END IF
914        END DO
915      END IF
916
917      IF (type_profil==4) THEN
918        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
919        DO j = 1, klon
920          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
921            rlatd(j))**4)
922          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
923            phy_sst(j, i) = 273.
924          END IF
925        END DO
926      END IF
927
928      IF (type_profil==5) THEN
929        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
930        DO j = 1, klon
931          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
932            *rlatd(j))**4)
933          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
934            phy_sst(j, i) = 273. + 2.
935          END IF
936
937        END DO
938      END IF
939
940      IF (type_profil==6) THEN
941        ! Méthode 6 "cst" valeur constante de SST
942        DO j = 1, klon
943          phy_sst(j, i) = 288.
944        END DO
945      END IF
946
947
948      IF (type_profil==7) THEN
949        ! Méthode 7 "cst" valeur constante de SST +2
950        DO j = 1, klon
951          phy_sst(j, i) = 288. + 2.
952        END DO
953      END IF
954
955      p = 0
956      IF (type_profil==8) THEN
957        ! Méthode 8 profil anomalies SST du modèle couplé AR4
958        DO j = 1, klon
959          IF (rlatd(j)==rlatd(j-1)) THEN
960            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
961              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
962            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
963              phy_sst(j, i) = 273. + x_anom_sst(pplus)
964            END IF
965          ELSE
966            p = p + 1
967            pplus = 73 - p
968            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
969              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
970            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
971              phy_sst(j, i) = 273. + x_anom_sst(pplus)
972            END IF
973            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
974          END IF
975        END DO
976      END IF
977
978      IF (type_profil==9) THEN
979        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
980        DO j = 1, klon
981          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
982            *rlatd(j))**4)
983          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
984            phy_sst(j, i) = 273. - 2.
985          END IF
986        END DO
987      END IF
988
989
990      IF (type_profil==10) THEN
991        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
992        DO j = 1, klon
993          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
994            *rlatd(j))**4)
995          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
996            phy_sst(j, i) = 273. + 4.
997          END IF
998        END DO
999      END IF
1000
1001      IF (type_profil==11) THEN
1002        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
1003        DO j = 1, klon
1004          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
1005            rlatd(j))**4)
1006          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1007            phy_sst(j, i) = 273.
1008          END IF
1009        END DO
1010      END IF
1011
1012      IF (type_profil==12) THEN
1013        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
1014        DO j = 1, klon
1015          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
1016            *rlatd(j))**4)
1017          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1018            phy_sst(j, i) = 273. + 4.
1019          END IF
1020        END DO
1021      END IF
1022
1023      IF (type_profil==13) THEN
1024        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
1025        DO j = 1, klon
1026          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
1027            rlatd(j))**4)
1028          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1029            phy_sst(j, i) = 273.
1030          END IF
1031        END DO
1032      END IF
1033
1034      IF (type_profil==14) THEN
1035        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
1036        DO j = 1, klon
1037          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
1038            *rlatd(j))**4)
1039          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1040            phy_sst(j, i) = 273.
1041          END IF
1042        END DO
1043      END IF
1044
1045      if (type_profil.EQ.20) then
1046      print*,'Profile SST 20'
1047!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1048
1049      do j=1,klon
1050        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
1051      enddo
1052      endif
1053
1054      if (type_profil.EQ.21) then
1055      print*,'Profile SST 21'
1056!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1057      do j=1,klon
1058        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
1059      enddo
1060      endif
1061
1062
1063
1064    END DO
1065
1066    ! IM beg : verif profil SST: phy_sst
1067    amn = min(phy_sst(1,1), 1000.)
1068    amx = max(phy_sst(1,1), -1000.)
1069    imn = 1
1070    kmn = 1
1071    imx = 1
1072    kmx = 1
1073    DO k = 1, year_len
1074      DO i = 2, nlon
1075        IF (phy_sst(i,k)<amn) THEN
1076          amn = phy_sst(i, k)
1077          imn = i
1078          kmn = k
1079        END IF
1080        IF (phy_sst(i,k)>amx) THEN
1081          amx = phy_sst(i, k)
1082          imx = i
1083          kmx = k
1084        END IF
1085      END DO
1086    END DO
1087
1088    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1089    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1090    ! IM end : verif profil SST: phy_sst
1091
1092    RETURN
1093  END SUBROUTINE profil_sst
1094
1095END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.