source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90 @ 3809

Last change on this file since 3809 was 3809, checked in by ymipsl, 10 years ago

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